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Some Remarks on Wild Observations* 


WituaM H. KrusKat** 
The University of Chicago 


Editor’s Note: At the 1959 meetings of the American Statistical Association held 
in Washington D. C., Messrs. F. J. Anscombe and C. Daniel presented papers on the 
detection and rejection of ‘outliers’, that is, observations thought to be maverick or 
unusual. These papers and their discussion will appear in the next issue of Techno- 
metrics. The following comments of Dr. Kruskal are another indication of the present 
interest of statisticians in this important problem. 


The purpose of these remarks is to set down some non-technical thoughts 
on apparently wild or outlying observations. These thoughts are by no means 
novel, but do not seem to have been gathered in one convenient place. 

1. Whatever use is or is not made of apparently wild observations in a statisti- 
cal analysis, it is very important to say something about such observations 
in any but the most summary report. At least a statement of how many observa- 
tions were excluded from the for. al analysis, and why, should be given. It is 
much better to state their values and to do alternative analyses using all or 
some of them. 

2. However, it is a dangerous oversimplification to discuss apparently wild 
observations in terms of inclusion in, or exclusion from, a more or less conven- 
tional formal analysis. An apparently wild (or otherwise anomalous) observation 
is a signal that says: ‘Here is something from which we may learn a lesson, 
perhaps of a kind not anticipated beforehand, and perhaps more important 
than the main object of the study.’’ Examples of such serendipity have been 
frequently discussed—one of the most popular is Fleming’s recognition of the 
virtue of penicillium. 

3. Suppose that an apparently wild observation is really known to have come 
from an anomalous (and perhaps infrequent) causal pattern. Should we include 
or exclude it in our formal statistics? Should we perhaps change the structure 
of our formal statistics? 

Much depends on what we are after and the nature of our material. For 
example, suppose that the observations are five determinations of the percent 
of chemical A in a mixture, and that one of the observations is badly out’ of 


* This work was sponsored by the Army, Navy and Air Force through the Joint Services 
Advisory Committee for Research Groups in Applied Mathematics and Statistics by Contract 
No. N6ori-02035. Reproduction in whole or in part is permitted for any purpose of the United 
States Government. 


** With generous suggestions from L. J. Savage, H. V. Roberts, K. A. Brownlee, and 
F. Mosteller. 
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line. A check of equipment shows that the out of line observation stemmed 
from an equipment miscalibration that was present only for the one observation. 

If: the magnitude of the miscalibration is known, we can probably correct 
for it; but suppose it is not known? If the goal of the experiment is only that 
of estimating the per cent of A in the mixture,-it would be very natural simply 
to omit the wild observation. If the goal of the experiment is mainly, or even 
partly, that of investigating the method of measuring the per cent of A (say in 
anticipation of setting up a routine procedure to be based on one measurement 
per batch), then it may be very important to keep the wild observation in. 
Clearly, in this latter instance, the wild observation tells us something about 
the frequency and magnitude of serious errors in the method. The kind of 
lesson mentioned in 2 above often refers to methods of sampling, measurement, 
and data reduction, instead of to the underlying physical phenomenon. 

The mode of formal analysis, with a known anomalous observation kept in, 
should often be different from a traditional means-and-standard deviations 
analysis, and it might well be divided into several parts. In the above very simple 
example, we might come out with at least two summaries: (1) the mean of the 
four good observations, perhaps with a + attached, as an estimate of the per 
cent of A in the particular batch of mixture at hand, and (2) a statement that 
serious calibration shifts are not unlikely and should be investigated further. 
In other situations, nonparametric methods might be useful. In still others, 
analyses that suppose the observations come from a mixture of two populations 
may be appropriate. 

The sort of distinction mentioned above has arisen in connection with military 
equipment. Suppose that 50 bombs are dropped at a target, that a few go wildly 
astray, that the fins of these wild bombs are observed to have come loose in 
flight, and that their wildness is unquestionably the result of loose fins. If we 
are concerned with the accuracy of the whole bombing system, we certainly 
should not forget these wild bombs. But if our interest is in the accuracy of the 
bombsight, the wild bombs are irrelevant. 

4. It may be useful to classify different degrees of knowledge about an ap- 
parently wild observation in the following way: 

a. We may know, even before an observation, that it is likely to be wild, 
or at any rate that it will be the consequence of a variant causal pattern. For 
example, we may see the bomb’s fins tear loose before it has fallen very far from 
the plane. Or we may know that a delicate measuring instrument has been jarred 
during its use. 

b. We may be able to know, after an observation is observed to be apparently 
outlying, that it was the result of a variant causal pattern. For example, we 
may check a laboratory notebook and see that some procedure was poorly 
carried out, or we may ask the bombardier whether he remembers a particular 
bomb’s wobbling badly in flight. The great danger here, of course, is that it is 
easy after the fact to bias one’s memory or approach, knowing that the observa- 
tion seemed wild. In complex measurement situations we may often find some- 
thing a bit out of line for almost any observation. 

c. There may be no evidence of a variant causal pattern aside from the observa- 
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tions themselves. This is perhaps the most difficult case, and the one that has 
given rise to various rules of thumb for rejecting observations. 

Like most empirical classifications, this one is not perfectly sharp. Some 
cases, for example, may lie between b and c. Nevertheless, I feel that it is a 
useful trichotomy. 

5. In case c above, I know of no satisfactory approaches. The classical approach 
is to create a test statistic, chosen so as to be sensitive to the kind of wildness 
envisaged, to generate its distribution under some sort of hypothesis of non- 
wildness, and then to ‘reject’ (or treat differently) an observation if the test 
statistic for it comes out improbably large under the hypothesis of nonwildness. 
A more detailed approach that has sometimes been used is to suppose that 
wildness is a consequence of some definite kind of statistical structure—usually 
a mixture of normal distributions—and to try to find a mode of analysis well 
articulated with this structure. 

My own practice in this sort of situation is to carry out an analysis both with 
and without the suspect observations. If the broad conclusions of the two analyses 
are quite different, I should view any conclusions from the experiment with 
very great caution. 

6. The following references form a selected brief list that can, I hope, lead 
the interested reader to most of the relevant literature. 


REFERENCES 


. C. I. Buss, W. G. Cocuran, ano J. W. Tuxsny, “A rejection criterion based upon the 
range,” Biometrika, 43 (1956), 418-22. 

. W. J. Drxon, ‘Analysis of extreme values,” Ann. Math. Stat., 21 (1950), 488-506. 

. W. J. Drxon, “Processing data for outliers,’’ Biometrics, 9 (1953), 74-89. 

. Franx E. Grosss, “Sample criteria for testing outlying observations,’’ Ann. Math. Stat., 
21 (1950), 27-58. 

. E. P. Kine, “On some procedures for the rejection of suspected data,” Jour. Amer. Stat. 
Assoc., 48 (1953), 531-3. 

. Jutius Lresiern, “Properties of certain statistics involving the closest pair in a sample 
of three observations,” Jour. of Research of the Nat. Bureau of Standards, 48 (1952), 255-68. 

. E. S. Pearson anp C. Cuanpra Sekar, “The efficiency of statistical tools and a criterion 
for the rejection of outlying observations,’’ Biometrika, 28 (1936), 308-320. 

. Paut R. River, “Criteria for rejection of observations,” Washington University Studies, 
New Series. Science and Technology, 8 (1933). 








TECHNOMETRICS Fesruary, 1960 


Statistical Estimation of the Gasoline Octane 
Number Requirement of New Model Automobiles 


CLAUDE S. BRINEGAR 
Rona.p R. MiLier* 
Union Oil Company of California 


A new method is proposed for estimating the octane number of gasoline required 
by new model automobiles. Tests of the assumptions underlying the method, and 
an illustration of its application, are given. 


SUMMARY 


An annual statistical problem in the automobile and petroleum industries 
is that of estimating, with confidence intervals, the octane number of the gasoline 
required to prevent knock in certain percentages (such as 50, 75, and 90) of 
the new model automobiles. The method currently in common use is to make 
nonparametric graphic estimates of the specified percentiles from smoothed 
curves of sample cumulative distribution functions. These functions are each 
based on about 25 observations of the octane number requirements of various 
major new model automobiles. Confidence intervals are computed from for- 
mulas appropriate for ranked estimates of percentiles of a normal distribution. 
In this paper the authors show that a simpler and statistically more efficient 
method is to estimate the required percentiles by computing z = ¢ + ks, 
with k selected from normal distribution tables. Confidence intervals can then 
be estimated from the approximate formula 


tw 1 as 
0, = 8 + om 


Examination of octane number data plotted on normal probability paper and 
detailed application of the x’ goodness-of-fit test provide strong evidence in 
support of the authors’ contention that the data conform closely to the normal 
distribution model. 


THE ENGINEERING PRoBLEM 
“Octane number” (abbreviated O. N.) describes the principal quality charac- 
teristic of a gasoline. In essence, the lower the number the more likely is it 
that the gasoline will give unsatisfactory performance in a modern automobile 
engine. The most common complaint is a distinct “knocking’’ during accelera- 
tion, an action that is both annoying to the driver and harmful to the engine. 
* Now employed by North American Aviation Company. 
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Technically, the O. N. of a given fuel is the percent by volume of iso-octane 
in an n-heptane-iso-octane blend that matches the antiknock quality of the 
fuel. (A special scale has been devised to measure numbers in excess of 100.) 
The O.-N. of the fuel in question is measured (“rated’’) by comparing it to 
known blends of n-heptane-iso-octane in standardized laboratory test engines 
under specified operating conditions. Though several methods are used, the 
most common one is the “Research” method. All octane numbers cited in 
this study are Research numbers. 

Each year, as the new automobiles are introduced, both the automobile and 
petroleum industries are interested in determining the O. N. of the gasoline 
needed to satisfy the requirements of the various new models. Surprisingly 
enough, although 100 new cars of model “X” may have the same factory- 
rated engines and transmissions, the O. N. of the gasoline needed to prevent 
knock in each of the 100 cars may (and often does) vary as much as 15 numbers. 
Thus, if refiner “Y’’ were marketing a 95 octane premium gasoline that would 
cause, say, 25% of all new model “X”’ cars to knock, the obvious result would 
be many dissatisfied owners of car ““X’’ and many unhappy customers of refiner 
“Y”. On the other hand, if refiner “Y’ marketed a 105 octane gasoline he 
would almost certainly meet the needs of all models of all new cars, but in order 
to do it he would be forced to raise the selling price in order to recover his higher 
manufacturing costs. In practice, most refiners strive for a balance of meeting 
the requirements of as many of the new cars as possible without unduly in- 
creasing their costs. A typical rule seems to be that of aiming at marketing a 
premium gasoline that “satisfies” at least 90% of each model of the major 
new cars. 

A non-profit corporation known as the Coordinating Research Council (abbre- 
viated CRC) each year undertakes the collection,: analysis, and publication 
of O. N. requirement data for the new model automobiles. The CRC (sponsored 
jointly by the Society of Automotive Engineers and the American Petroleum 
Institute) obtains, from about 30 cooperating oil company and automobile 
research laboratories, several hundred reports on the O. N. requirements of 
the major new cars. Considerable pains are taken to see that each field report 
is prepared under the same controlled conditions. In this way it is hoped that 
the data on each new model car may be regarded as random observations from 
a homogeneous population. Since it costs several hundred dollars to satisfac- 
torily estimate the O. N. requirement of a single car, the need for cooperative 
action is evident. 

The authors became interested in the statistical treatment of the CRC data 
when asked certain questions about it by members of Union Oil’s research 
department and when the statement was noted, in the CRC reports for 1955 
and 1956, that “‘the exact statistical method of treating this type of data has 
not yet been decided upon ....’” [4, Appendix C]. 


THE STATISTICAL PROBLEM 


Briefly, the present CRC method of analysis is as follows: Suppose 23 reports 
on car “X” (all with the same type engine and transmission) are submitted 





re woe ewe eh 


OS = § ODO mee me ee EG hl 


STATISTICAL ESTIMATION OF GASOLINE OCTANE 7 


by the cooperating laboratories. The raw observations are as shown in Table 1. 

The data are first ranked, and the percentiles computed by dividing the 
ranked order by n + 1. This is shown in Table 2. (The form i/n + 1 is used 
to provide unbiased estimates of the areas above and below points on the cumu- 
lative distribution of octane numbers. See Mood [9, pp. 385-88].) The per- 
centiles are then plotted and a smooth curve is fitted free-hand, with some 
attempt generally made to go through the center of the “clusters.” Chart 1 
shows one such curve for model ‘‘X”’. (The CRC reports designate the “per- 
centiles” as “percent satisfied” and plots them on the abscissa instead of the 
ordinate; here we follow the more traditional approach.) For purposes of the 
CRC study, the O. N. requirements corresponding to the 50th and 90th per- 
centiles are read directly from the curve. One of the most important parts of 
the CRC report is a table giving these 50th and 90th percentile estimates for 
each of the 15 to 20 models included in the annual survey. In addition, the 
charts showing the actual smoothed curves are also published. From Chart 1 
we can read the 50th percentile as 90.0 and the 90th as approximately 93.5. 


TABLE 1 
Vehicle Code—“‘X”’ 


O. N. Requirement 
Sample Number (Primary Reference Fuel)* 


1 
2 
3 
4 
5 
6 
7 
8 
9 


* Two basic fuels (“Primary Reference” and “Full Boiling Range’’) are used in the CRC 
studies. Throughout this study we limit ourselves to the Primary Reference fuel. Results 
obtained from the two fuels have been found to be very closely correlated. 
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TABLE 2 
Ranked Order and Percentiles 


O. N. Percentiles 
Ranked Order Requirement (“% Cars Satisfied’’) 


82 
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87 
91 
95 


This means that a 90.0 octane gasoline will prevent knock in about 50% of 
the new model “‘X”’ cars and a 93.5 octane gasoline will do the same for some 90%. 

The CRC reports relegate the problem of confidence intervals to an appendix. 
Appendix C of the 1957 report notes that for purposes of computing confidence 
intervals, “it was assumed that the car samples were randomly selected and 
that the data obtained were normally distributed” (p. 127). The formulas used 
for the standard deviations of the 50th and 90th percentiles are the following 
(see, for example, Kendall [7, p. 238]): 


Op, = 1.25330/+/n 
Oy.e.0 = 1.70940/+/n. 


Though not noted in the CRC reports, there is the implicit assumption that 
the graphic estimates of the percentiles are the same as the unbiased estimates 
computed directly from the ranked data. 

After analysis of reports and recommendations derived from the CRC O. N. 
data and discussion with research chemists, engineers, and others who use the 
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OCTANE NUMBER REQUIREMENT.CURVE, 
MODEL "x", 1955* 
100 


PERCENT 


80 


88 92 
OCTANE NUMBER REQUIREMENT 


#% COMPUTED BY D'ivIDING RANKED ORDER By n+ 


data, the authors believe that the CRC method of making non-parametric, 
graphic estimates of the octane numbers is not the best available one. We 
believe that it can be demonstrated that the O. N. data are, in fact, essentially 
normally distributed and that, consequently, a better method is the following: 
Compute the sample mean (#) and standard deviation (s), and then estimate 
the required percentiles by computing z, = ~ + ks, with k selected from normal 
distribution tables. Commonly used values of k are: 


Percentile k 


5 —1.645 
10 —1.282 
25 —0.674 
50 0 
75 +0.674 
90 +1.282 
95 +1.645 


Confidence intervals are easily computed from the approximate formula for 
the standard deviation of z [6, p. 58]: 


1 k? 
o,, = Sa/—- 


a’ 1) 


Assuming that normality (or near-normality) can be established, this method 
has several advantages, in the authors’ view, over the present CRC method. 
These include better statistical efficiency (for example, compare the above 
variances for the 90th percentile and for + 1.282s); ease of calculation, under- 
standing, and explanation; ease of adding later observations to partial data 
without complete recalculation; and reproducibility. 
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TaBie 3 
Comparison of Three Estimates of 50th and 90th Percentiles of O. N. Requirements* 
A. 50th Percentile Estimates 
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B. 90th Percentile Estimates 


1956 
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97.4 97.6 
96.5 96.9 
100.0 100.3 
97.0 96.7 
93.0 92.4 
99.0 98.8 
97.5 98.8 
96.0 96.0 
96.0 $5.9 
91.6 91.1 
97.5 97.4 
94.0 94.0 
92.0 92.3 
96.5 96.1 
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* “P,” designates graphic estimates made by CRC; “P,’, estimates made by interpolation 
of ranked data; “‘z’”’, estimates made by computing # and # + 1.282(s). The codes “D”, “G”, 
etc., refer to CRC model designations. 
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We can illustrate its use by computing the 50th and 90th percentiles for the 
23 observations for ““X’”’ shown in Table 1: 


n 3 
n se 
2 = i=1 aa 2 92 
nn — ), 7 


50th percentile: 89.91 + 0(2.92) = 89.91 


90th percentile: 89.91 + 1.282(2.92) = 93.65 


Two-sigma confidence intervals are 89.91 + 1.22 and 93.65 + 1.66. (This 
compares with the corresponding two-sigma confidence intervals of 90.0 + 1.53 
and 93.5 +°2.08 computed by the present CRC method.) 

The entire cumulative distribution curve can be easily plotted. on either 
traditional coordinate paper or, perhaps most simply, as a straight line on 
arithmetic probability paper. The confidence interval band also can be directly 
plotted. 

In the following section we investigate the basic question of just how satis- . 
factory is the assumption of normality. 


Test or Hypotuesis THAT O.N. DATA ARE NORMALLY DISTRIBUTED 


Before presenting evidence in support of the normality assumption, it is of 
interest to look briefly at the kinds of estimates obtained (1) by the graphic 
method used by the CRC (described previously), (2) by estimating the per- 
centiles directly from the ranked data, and (3) by using z, = # + ks, the method 
proposed by the authors. Table 3 summarizes these three estimates of the 50th 
and 90th percentiles for the 1955-57 octane data. 

It is clear that in a general way the three estimates agree quite well. The 
averages are in close agreement and individual variations of as much as one 
O. N. are unusual. One interesting conclusion comes from study of the relation- 
ship between the graphic and ranked estimates of the two percentiles. It appears 
that the curves fitted freehand to the data have a tendency to be influenced 
too much by the upper extreme observations in the sample. This can be seen 
by comparing the signs of the differences between “‘P,’”’ (graphic estimate) 
and “P,” (ranked estimate). 


50th Percentile 90th Percentile 
Pi > P, 13 25 
P, = P; 23 8 
P, <P, 8 11 


Thus, while the split at the 50th percentile is a reasonable random fluctuation, 
the split at the 90th is not (see, for example, Dixon and Massey [5, p. 417]). 
Two separate methods have been used to investigate normality of the basic 
O. N. data. These are (1) studying the behavior of the individual cumulative 
distribution curves when plotted on normal probability paper and (2) testing 
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the ‘“‘goodness of fit’’ of each set of data to the normal hypothesis by use of 
the chi-square (x’) test. We will discuss the results in turn. 


NorMAL PRoBABILITY PAPER 


While there is no simple statistical method available for testing the hypothesis 
that a sample plotted on normal probability paper is or is not a random varia- 
tion from a straight line (the normal distribution hypothesis), examination of 
the plots of several independent samples should indicate if any obvious con- 
sistent non-linear tendencies exist. Chart 2 shows such a plot for the 1956 CRC 
data. There we have arranged the 14 models included in the 1956 survey so 
that none of the cumulative distribution curves overlap. (The extreme tails 
are omitted since their wider fluctuations detract from our basic purpose of 
studying linearity in the range of main interest.) Inspection of this chart, in 
the authors’ view, reveals no obvious consistent departures from linearity. 


CHART 3 


AVERAGE OF 
OCTANE NUMBER REQUIREMENT CURVES 


COMPUTED BY TWO METHODS, CARS IN 1955-57 SURVEYS 


PERCENT 


w 
2 
@ 
2 
” 
= 
2 
o 
a 
o 
° 
= 
a 
2 
w 
+ 
x 


( ari 


‘ 
1956 1957 


————_ Computed by Norma! Distribution Method 
———— Computed trom Curves such as shown in Chart 2 


OCTANE NUMBER REQUIREMENT 
(VARYING SCALE} 


Chart 3 shows the result (as a dotted line) of averaging horizontally the 14 
curves in Chart 2, and also the averages for the 13 curves of 1955 and the 17 
curves of 1957 (computed but not shown individually). The solid straight lines 
are the averages of the individual percentile estimates computed directly by 
the normal distribution method. Comparison of the solid and dotted lines 
suggests one non-linear bias: The slight tendency for the lower tail-of the empir- 
ical distribution curve to rise a little more rapidly than the normal distribution 
model. From the 50th percentile on, however—the part of the curve of greatest 
interest—the dotted and straight lines agree quite well. 

Thus, while the behavior of the O. N. data on normal probability paper 
does not establish beyond doubt that the normal distribution model is a suitable 
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one, it strongly suggests (by “weight of evidence”) that if the underlying distri- 
bution does depart from normality, the departure is slight, probably confined 
to the lower tail. 


Cui-SquaRE Test of GOODNESS OF Fit 


Using x’ to test normality of small-sample data is, unfortunately, borderline 
between an art and a science. This is clearly brought out by W. G. Cochran 
in two detailed articles on the elementary uses of the x’ test [2, 3]. Cochran 
points out that while rules commonly appearing in elementary statistics books 
recommend a moderate number of equal-width class intervals (say between 
10 and 25) and the grouping of the “tails” so that the minimum expected value 
is five or ten, the “recommendations of those who have looked into this problem 
are contrary to current practice” [2, p. 332]. Mann and Wald [8] recommend 
selecting class intervals so that the expected numbers are about the same, and 
prescribe a method for selecting the optimum number of classes for large 
sample sizes. Williams [10], in an elementary description of the Mann and 
Wald method, notes that it ‘may be valid for sample sizes as low as 200 and 


TABLE 4 


Example of x? ‘‘Goodness of Fit’’ Test 
Model “X”, 1955 


a. Five Class Intervals 


Class Interval Expected Observed 


Up to 87.5 4.68 
87 .5-89.5 3.17 
89.5-90.5 8.77 
90.5-92.5 2.03 
92.5 and over 4.35 


Total 23.00 


b. Divide First Class Interval 


Up to 85.5 1.51 2 
85.5-87.5 3.17 1 


es a 


Total 23.00 23 


c. Divide Fifth Class Interval 


92.5-94.5 3.02 
94.5 and over 1.33 


Total 23.00 


d. Seven Class Intervals 
Class intervals 1 and_2 of b + 3, 4, and 5 of a + 6 and 7 ofc 
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may be true for considerably smaller samples” (p. 83). Though Williams does 
not hazard a guess as to the absolute minimum, the authors have inferred that 
it is greater than the samples in this study. 

The O. N. data for the years 1955-57 are in 44 separate sets, ranging from a 
high of 37 observations to a low of 16, with an average sample size of 27. Since 
each set has presumably different (and unknown) means and variances, the 
raw data cannot be combined. However, once the x’ tests have been made, 
the separate x” values can be pooled to provide a single test of the normality 
of all sets taken together. This is the procedure we follow. 

The class intervals pose something of a problem. If we select them so that 
each expected value is not less than about five, it is clear that five intervals 
are the absolute minimum. Cochran [3, p. 420], however, believes that in testing 
the fit of unimodal distributions (such as the normal or Poisson), it is quite 
satisfactory to group the tails so that the minimum expectation at each tail 
would be not less than one. Such a grouping would obviously provide a more 
sensitive test of the tails for our small samples than one limited to an expected 
value of five. 

After considering the various methods of applying the x’ test to the O. N. 
data, the authors decided upon the following course of action: 

1. Initially, all sets would be separately and arbitrarily grouped into five 
class intervals. The first and fifth intervals would be wide enough so that each 
would include about one-fifth of the observations. 

2. Second, the first and fifth intervals would be each split into two separate 
intervals, so the smallest expectations would be not less than one. 

3. The x’ tests would be applied to the four grouping possibilities: 


a. Initial intervals one, two, three, four, and five. 

b. Split interval one, and initial intervals two, three, four, and five. 

c. Initial intervals one, two, three, four, and split interval five. 

d. Split interval one, initial intervals two, three, four, and split interval 
five. 


Since the degrees of freedom for the x’ test for normality are equal to the number 
of class intervals less three [5, p. 227], the degrees of freedom for these variations 
are, respectively, two, three, three, and four. It was believed that these alterna- 
tives would (1) provide some indication of the effect of small expectations on 
the x’ test and (2) indicate whether the fit of the lower tail was poorer than 
that of the upper, as Chart 3 suggests. Table 4 illustrates these calculations 
for the set of data shown in Table 1. 

Table 5 shows the average expected value of each class interval for the four 
possible groupings. 

Table 6 summarizes the actual x’ tests for the 44 sets for the four combinations. 
As shown at the bottom of this table, none of the pooled values of x’ is signifi- 
cant at the 5% level. Thus, we have additional support for the conclusion reached 
on the basis of our analysis with normal probability paper that the O. N. data 
may be accepted as normally distributed. Since the “seven interval” test provided 
the closest fit to the normal model (probability of about 0.3 of a worse fit), 
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TABLE 5 
Average of Class Interval Expected Values 


Six Class Intervals 


Five Lower Tail Upper Tail Seven 
Intervals Divided Divided Intervals 


ee , 


Total 27.5 : 27.5 


we have further evidence of the validity of Cochran’s statement that expected 
values of about one in the “tails” are satisfactory in testing for normality. 
Though difficult to judge statistically, we believe that it is also significant that 
the six-interval test with the split upper tail produced a better fit than the 
one with the split lower tail (probability of “worse” fit of about 0.25 as com- 
pared with about 0.10). This supports the conclusion reached by inspection 
of Chart 3 that the departure from normality is mainly confined to the lower 
tail—the one of least interest insofar as the O. N. data are concerned. 

Three final comments are in order. First, as Table 6 shows, model L-2 for 
1957 gave far and away the poorest fit of the 44 sets. Inspection of the individual 
data for L-2 suggests that these observations are actually from two overlapping 
normal distributions (the large value of x’ results from the absence of observa- 
tions near the sample mean). Field reports have noted that this model is particu- 
larly susceptible to “rumble”, a new form of knock associated with high-com- 
pression engines. It is possible that because of the difficulty of objectively 
measuring rumble, the 28 observations on L-2 might have considerable non- 
random bias. Conceivably, by further digging we might produce sufficient 
justification for excluding this “maverick”. Such action, of course, would sub- 
stantially improve the x’ fit of the remaining 43 sets. 

Second, if the apparent slight departure from normality in the lower tail 
is real, then it is probable that by accumulating more and more observations 
we should eventually reach a significant value of x’. This difficulty (i.e., the 
cumulative effects of small differences between the hypothesized and actual 
distributions) has been discussed in the literature, though no general conclusions 
seem to have been reached. A colorful comment on the problem was made by 
Berkson [1]: “For we may assume that it is practically certain that any series 
of observations does not actually follow a normal curve with absolute exactitude 
in all respects, and no matter how small the discrepancy between the normal 
curve and the true curve of observations, the chi-square P will be small if the 
sample has a sufficiently large number of observations in it’ (pp. 526-27). 

The third comment concerns this question: If the O. N. data are approximately 
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TABLE 6 
Summary of x? Tests for O. N. Data 1955-57 


Six Class Intervals 


Five Class Lower Tail Upper Tail Seven Class 
Year and Code Intervals Divided Divided Intervals 


1955-D ; , é 5.25 
G 
H 
I 
L-2 
M 
N-A 
N-M 
O-1 
O-M 


P 


Total 


Degrees of 
Freedom 
P(x? > Total) 


2. 
a 
Z. 
6. 
4. 
4. 
2. 
3. 
¥. 
2. 
0. 
6. 
0. 
0. 
3. 
3. 
3. 
E. 
E; 
2. 
E; 
2. 
4, 
2. 
EB: 
6. 
3. 
2. 
2. 
2. 
12. 
2. 
2. 
3. 
4. 
2. 
3. 
2. 
6. 
40. 


19 
59 
78 
a 
75 
54 
40 
57 
88 


es 


_ 
(JS) 
bo 
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normally distributed, as the authors believe they have shown, why has this 
not been recognized before? In early discussions of this statistical problem 
with several chemists and engineers who work closely with O. N. data, the 
authors were told repeatedly that the basic data could not be normally distri- 
buted because octane numbers are not of “‘equal size” as we move up the scale. 
That is to say, it takes only a slight modification in engine adjustment or design 
to boost the requirement several numbers in the lower range, whereas a large 
change is needed in the upper range. Granted the technical accuracy of this 
statement, how then do we account for the fact that as O. N. requirements 
have increased with time (especially between 1955 and 1956, as can be seen 
in Table 3) the normal distribution has continued to be a satisfactory model? 
The answer may lie, at least partly, in the fact that there is not just one, but 
rather several determinants of a given car’s O. N. requirement. Thus, if we 
start with 100 brand new, carefully assembled cars with identically-designed 
engines and transmissions, we might find that the distribution of O. N. require- 
ments would have a non-normal shape, with somewhat smaller dispersion than 
we have found in the data we have studied. However, to this basic design re- 
quirement we must add the effects of the engine’s factory adjustment and 
subsequent tune-ups, the owner’s driving habits, and the type of gasoline and 
lubricating oil that has been used. This action could result in a total cumulative 
effect that would permit the central limit theorem to be applied, thus creating 
a reasonably normal composite distribution (see, for example, [9, pp. 136-39]). 
While the authors cannot present data in support of this hypothesis, it does 
at least suggest the possibility of a theoretically valid explanation of what we 
have, in fact, observed. 
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The paper examines what can be said about overall grand lot quality based upon 
random samples from batches or sublots disposed of sequentially as they are produced. 
Assurance always depends upon how the sample is selected and upon assumptions 
(restrictive or otherwise) about the population from which the sample came. The 
wisdom of assuming that a statistically controlled process is a single process is ques- 
tioned. Moreover, the paper points out properties of the sample, stratified as it is 
over time of production, and of the population it represents, decimated as it may 
be by rejection of some sublots, that make it unwise to expect certain assurance 
statements from a vendor based on such samples. The paper demonstrates a valid 


method for making what quality assurance statements are possible based on such 
samples. 























1. INTRODUCTION 


Suppose a vendor is asked to produce a large number of fuel elements for 
a reactor. In order that production of subassemblies of fuel elements not be 
held up, individual lots of elements are released or rejected as they are produced, 
at least partly on the basis of evidence contained in destructive measurements 
on a simple random sample of n elements from each lot of size N. When, however, 
the total required number of elements has been accepted, what assurance is 
provided by the stratified simple random sample of destructively measured 
elements about the quality of accepted elements? 

Fuel elements are expensive. Consequently n, the number of sample elements 
destroyed per lot, is kept small, possibly as low as 1. On the other hand, every 
effort is made to maintain strict control over a production process that was 
proven during pre-production qualification to be capable of producing a high 
proportion of acceptable fuel elements. However, a =rocess continuously in 
statistical control throughout total production may not be a process with 
constant parameters, a point particularly relevant to the problem of estimating 
the quality of total production and to setting confidence limits on that quality. 
To estimate and to establish confidence limits for the quality of total accepted 
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product, we make no assumptions whatever about the nature of the process. 
We will not even assume that a single process produced all accepted lots. Such 
assumptions add nothing to our ability to set confidence limits on the quality 
of the elements contained in our specific accepted lots each of size N — n. 

To make our problem precise, suppose that K sublots each of size N are to 
be produced and that K(L) sublots (or, more correctly, sublots diminished by 
the number of sampled elements) constitute a reactor load [K(Z) < Kj]. Of 
these K sublots K, are accepted and K, are rejected; K, and K, are, of course, 
not known in advance. If K, > K(L) no further sublots are involved. If K, < 
K(L), further sublots are produced until K(Z) — K, additional sublots are 
accepted. The problem is to evaluate the quality of the K(L) accepted sublots 
each of size N — n. 

We have intentionally chosen to speak of quality rather than percent de- 
fective among accepted elements, for certain peculiarities arise in consideration 
of the inferences that can be made from our stratified random sample as opposed 
to a simple random sample from accepted lots. Consider the K lots sure to be 
produced. If our sample of nK destructively measured elements were a simple 
random sample from the NK elements, the sample order statistics could be used 
to estimate the percentiles of the total production, accepted and rejected, and 
to get confidence limits on those percentiles. 

Our sampling is not simple random, but stratified random. In view of this, 
what can be said about percentiles of the distribution of elements in all lots 
produced? It is possible to estimate the percentile of each lot separately. But, 
an ordering of the total production may in no sense coincide with the ordering 
of lots. Valid confidence limits on a percentile of total production could be 
obtained by getting confidence limits on the same percentile of each lot and 
taking the outermost limits for total production. This, however, could be practi- 
cally useless since the sublot confidence would have to be taken at (1 — a’) 
where (1 — a’)* = 1 — a, in order for the confidence on total production to 
be 1 — a. The fact that n is quite small only makes matters worse. This appears 
to be the best we can do by way of assurance from variables information when 
we do not know that one and only one process produced all lots. If only accepted 
lots are considered there is a further complication: if acceptance of sublots 
is based on sample evidence of defectives, the evidence from accepted lots is 
from censored samples. Matters get worse instead of better, when an additional 
K(L) — K, accepted lots must be considered (see Section 4). Thus, we cannot 
as a practical matter (although perhaps it is feasible in principle) make even 
the weak statement indicated above for percentiles of the distribution of accepted 
elements. 

Variables information can be compared to some previously assigned upper 
or lower limit, or attributes can be used directly to decide whether a destructively 
measured element is acceptable or not. Then, to a good approximation, we can 
get confidence limits for the percent defective among accepted elements, either 
for K, lots from K lots, or for K, out of K lots plus K(L) — K, additional lots 
in case K, < K(L), completely independently of the quality of individual 
lots. These limits will tend to be practically useful only when (K(L) — K,)/K 
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and K,/K are small; that is, for generally good processes or possibly combina- 
tions of generally good processes. 


2. SUMMARY 
The following notation is used: 


P, = true proportion of defectives among all accepted lots. Each lot is 
diminished by its sample size 
K(L) = the number of lots (each diminished by its sample size) necessary 
to give a required total number of items. 
K = the number of lots sure to be produced to fill the given order [specified 
in advance and > K(L).]. 
K, = the number of lots accepted out of K. 
nm; = sample size from ‘‘2’’th lot. 
r; = sample number of defectives in ‘’’th lot. 
N, = number of elements in “‘2’’th lot. 
i, = standardized normal deviate exceeded with (one sided) probability a. 
f = sampling rate per lot (0 < f < 1). 


For simplicity, the K, accepted lots out of K are designated as the Ist, 
2nd, --- , K,th lots. If K(Z) > K, , the subsequently produced and accepted 
K(L) — K, lots are designated as the (K + 1)th, (K + 2)th, --- ,(K + K(Z) — 
K,]th lots. The rejected lots out of the first K lots are designated as the 
(K. + 1)th, --- , Kth lots. We also define 3 


Ji+ea-9 / EN. | 


This paper shows that an approximate and conservative 100 (1 — a)% upper 
confidence limit for P, is given by 


2.1) Q= 


K+K(L)-—Ka K+K(L)-Ka 


Ka 
Q+ pe Ne a oo x ri 
(2.2) ee eee 
ee p| Sm, +. 2 w, | 

There are two approximations involved in (2.2). The important one is that of 
approximating the distribution of the observed total number of defectives in a 
reasonably large number of independent random samples from a series of lots 
by the normal distribution. The second approximation is the trivial one of 
replacing a factor 1/(N; — 1) in the variance of the observed number of defectives 
in the “#’th lot by 1/N; . The conservatism of (2.2) stems from replacement 
of the variance of observed total number of defectives by a somewhat larger 
number; this can be shown (non-asymptotically) to increase the confidence 
coefficient to be ascribed to (2.2). 
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3. AN EXAMPLE 


Assume an order with a vendor for approximately 225 (K) lots of 160 (N) 
elements each. The vendor is required to destructively sample a minimum of 
one element per lot but is also required to furnish at least 95% assurance that 
no more than 1% of the total of accepted elements are defective. 

A simple random sample of 225 elements containing no defectives from the 
grand lot of 36000 thoroughly mixed elements would yield about, but not quite, 
the required assurance. 

To insure safe reactor performance only fuel element manufacturing processes 
with extremely low proportion defectives are ever released for production. In 
order that we might expect a few defectives, suppose, strictly by way of example, 
that the process as qualified produced 3% defectives. Further suppose that 
following pre-production and for the first 3 of lots produced, the probability 
of a defective is }% and for the remaining 3 of lots, the probability of a defective 
is $%. We would then expect 60 defectives or 3% among the grand lot of 36000 
elements. 

Of the first 75 lots submitted suppose 

No. of lots No. of defectives per lot 


50 0 
20 1 
4 2 
1 3 


Of the last 150 lots submitted suppose 
No. of lots No. of defectives per lot 
123 0 


25 1 
2 2 


If the vendor samples at the minimum rate of 1 (n) per lot 
No. of the 20 lots with 1 def. Prob. of getting the def. in the sample 

0 - 8822 

i .1110 

2 .0066 

3 .0002 


No. of the 26 lots with 1 def. Prob. of getting the def. in the sample 


0 -8550 
i . 1344 
2 .0101 
3 .0005 


No. of the 4 lots with 2 defs. Prob. of getting a def. in the sample 


0 : -9510 
1 0481 
2 .0009 


No. of the 2 lots with 2 defs. Prob. of getting a def. in the sample 


0 9752 
1 .0247 
2 .0001 
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No. of the lot with 3 defs. Prob. of getting a def. in the sample 

0 -9813 

1 .0187 
The average number of lots from which we expect to get a defective in the sample 
is only .37. So even with our proportionately high process proportion defectives 
not a single defective might appear in any of the 225 samples. Taking r; = 0, 
+= 1,2,---,Kand K, = K, (2.2) gives a confidence limit for the proportion 
defective in accepted lots: 


£ ’ -1 
225 (1 ~ &) ; 
For a = .05, an upper confidence limit for the proportion defective in accepted 
lots is P, < .012, with confidence at least .95. 


4. DERIVATION OF RESULT 


Let P, be the proportion of defective elements in the first lot before the sample 
is taken; let P, be the proportion defective in the second lot before sampling 
and in general let P; be the proportion defective in the “‘z’’th lot before sampling. 

In the ‘2’’th lot, P; is estimated by r,/n where r; is the number of defectives 
in our sample from that lot. 

If only the K lots sure to be produced are considered, then each r; inde- 
pendently has a hypergeometric distribution with mean NP; and variance 
[n(N — n)/(N — 1)]P.(1 — P,). Consequently, if K is large or if n is large 
but < N, or both, we can assert that, approximately, the quantity 


Kx x 
YSr—n > P; 
1 1 


Ve =n Pil — P) 


is normally distributed with zero mean and unit standard deviation. 
To this approximation 


(4.2) PriZ > —t,} = 1—a, 


where ¢, is the standardized normal deviate exceeded with probability a. 
However, from the sure inequality 


K x 2 
(4.3) LP - (> P.) /K > 0, 
1 1 
it can be shown to the same approximation that the quantity 
K K 
> <= - P; 
1 1 
N K K 2 
vibe. - (Lr) «| 


1 1 


(4.1) Z= 


(4.4) Zt = 


is such that 
(4.5) Pr{Z* > —t.} >l—a. 
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Solving this last inequality for >>* P; , it follows that, since Z* is monotone 
decreasing in )-* P; , 


wit; N—n te , te van | No —n te 1(> y] 
on Sp, i Do ety tte Pt Dre dr 
N—n ta 


+N aK 


with confidence at least 1 — a. For simplicity only, suppose the sublots accepted 
were the first K, . Then the true proportion defective in accepted sublots (after 
sampling) is 


(4.7) -, = z bP ee 3 run a — n/N). 
Denoting the right hand side of (4.6) by U 
(4.8) PS Llu - = P,— Yr, x a — n/N). 


The right hand side of (4.8) will only be increased if we take P;(¢ = K, + 
1, --- , K) equal to zero. Hence, for K, accepted sublots out of a certain pro- 
duction of K sublots 


(4.9) PS Llu - ¥ r un \a — n/N) 


with confidence at least 1 — a. 

Now suppose K, < K(L) so that additional sublots are produced. Denote 
the proportion defective of the additional accepted lots before sampling by 
Pri, +++ , Pesxct)-x, and the observed defectives in samples of n from each 
accepted lot by resi , Tx+2 » *** » 1xK+K(L)-Ke « Then an upper confidence limit 
is derived from the confidence statement on )_* P; by writing 


K+K(L)—Ke K K+K(L)—-Ka 


Ka 
pee ee i ee 


K+1 < atl K+1 

(4.10) K(L) K(L) 
1 a 

so that for the proportion defective in accepted lots after sampling, P, , 


P, U + {KL — K.} — Sinn 


" < EH l 


(4.11) 


K+K(L)-Ke 


dN a - 9/m) 


One might ask why we did not obtain limits in this latter case by considering 
the total number of lots, K’ say, that it was necessary to produce in order to 
obtain K(L) acceptable lots in the same fashion as we treated the K lots which 
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we were sure would be produced. The reason is that K’ is a random variable so 
that the theory described above would not hold; it might be added that a treat- 
ment taking the random character of K’ into account leads to quite intractable 
theory. Thus, as suggested earlier, if only a1 good process is qualified and if 
only moderate departures therefrom occur, K — K, and K(L) — K, should 
be small and therefore the confidence limit will be useful. However, no matter 
how large K — K, or K(L) — K, may be, the confidence limit is valid to the 
extent of our normal approximation. 

For K at all large and n < N the confidence limit can be simplified by replacing 
U by 


K r: £ t vé K 1 K 2 
Pa, =e aa == —# —_—— ; 
(4.12) Ut = ee tet at Dn (dx). 


It is easy to generalize the above result, to essentially the same degree of 
approximation, to the situation where the sublots are not all the same size but 
proportionate sampling is done within lots. That is, a constant proportion, f, 
of elements is destructively sampled from each lot. Let P; be proportion defective 
in the “z’’th sublot before sampling, NV; the lot size and n; the sample size and 
assume, for simplicity, that only the first K, sublots were accepted. Then the 
true proportion defective in accepted sublots (after sampling) is given by 


Ke 


Ka 
> NP; baa dr; 


(4.13) P, = 


(4.14) 


is distributed approximately unit normally under the same conditions as already 
noted. Arguing essentially as before to the degree of our approximation, 


5, Se 
SSS 2 
Tee - (eae) e] 


Confidence limits on >>** N;P; follow from this result in an obvious way. 
Although the discussion above has been in terms of destructive sampling, 
the same type of argument is applicable to non-destructive sampling with or 
without rectification of sample defectives. The only alteration necessary would 
be in the definition of the true proportion of defectives in accepted lots. 
One further point is worth noting. After any number of lots < K have been 
produced one might wish to know whether the evidence already at hand precludes 


(4.15) fs 
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the possibility of obtaining the required assurance. This could easily be deter- 
mined by assuming that future sublot samples would show no defectives and, 
if we were appraising what might happen before K sublots had been produced, 
that at least all future sublots up to the Kth were accepted. Use of (4.11) or 
(2.2) would then yield the most optimistic estimate of the assurance that could 
be expected. Of course, one would prefer to be able to decide at any point in 
time before completion of production that it is certain that the desired assurance 
will be obtained; unfortunately, and obviously, this is impossible. 
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Elements of the Theory of Extreme Values’ 
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Wayne State University and Stanford University 


The theory of extreme values finds many applications in modern technology. 
This paper is an attempt to unfold the basic elements of the theory. Frequent ex- 
amples are given to illustrate both derivations and application. 


1. INTRODUCTION AND SUMMARY. 


The statistical theory of extreme values has, in recent years, been very useful 
in connection with a variety of applied problems (see references). However, 
the engineer faced with the task of learning the elements of this theory and 
then applying it, finds that the literature is both extensive and difficult. 

While the theory of extreme values, viewed in its entirety, is a beautiful and 
in places esoteric theory, its basic elements can be given in a fairly compact 
form. An attempt at such a presentation is given in this paper. In particular 
we show how various important limiting distributions can be obtained. We 
have adopted the device of illustrating various theoretical points by means 
of appropriate examples. Toward the end of the paper, some applications of 
extreme value theory are mentioned. 

A detailed theoretical treatment of the subject together with indication of 
applications of the theory is given in Gumbel [7]. 

Throughout this paper we center attention on the following problem: 

Consider samples of size n, drawn independently and at random from a 
population described by its probability density function f(x) or alternatively 
by its cumulative distribution function, F(x) = f%.. f(y) dy. Let the observations 
drawn be {2z, , 22 --- , Z,}. Define the two statistics: 


Yn = min {x, , 22, °** Ln} 
and 


2, = max {2 , 2%, °*** » Zn}. 


The problem is to find the density and cumulative distribution functions of y, 


and z, , respectively, and to give suitable asymptotic formulae for these functions 
when n gets large. 


2. DisTRIBUTION OF SMALLEST VALUES. 
The distribution of y, is found as follows: 
(1) Prob (y, > y) = Prob (all z; > y) 
= Prob (1, > y, 22 > y,°°',%. > y) = [1 — FY. 
1 Research sponsored in part by the Office of Naval Research. 
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Therefore, 


(2) G,(y) = Prob (y, < y) = 1 — Prob (y, > y) 
=1— [1 — Fy). 


This is the cumulative distribution function of y, . The density function g,(y) 
is found by differentiating G,(y). This gives us 


(3) gly) = Gry) = nfty[l — Fywr". 


Let us consider some special cases in which we get a simple result for any value 
of n, the sample size. 


Example 1: Underlying distribution is uniform over [0, A]. In this case 


(4) 


fa)=4, O<2<A 


= 0, elsewhere 


F(z) t<0 
e6s24<. A 


= TPs A, 


In this special case G,(y) becomes 


(6) G.y)=0, y<0O 


jue tit ac 
= 1 ae ? O<y<A 


Ly . are 


(7) gy) = Gily) =n Lo 
= 0, elsewhere. 
Example 2: Underlying distribution is exponential. In this case 
(8) f(x) = re, x20 


= 0, elsewhere 


and 
(9) . «<<, 
1-—e™“, 2#>0. 
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In this special case 
(10) G.y)=1-e™, y2O0 


= 0, elsewhere 
and 


(11) gy) = Gity) =mre™, y>O0 
= 0, elsewhere. 
Note that if x is distributed exponentially with parameter \, then 
Yn = Min (4%, , X2, °** 5 Ln) 


is distributed exponentially with parameter ni. 

While the G,(y) and g,(y) in Examples 1 and 2 were of a particularly simple 
form, these functions are in general very complicated. For example, if the 
underlying p.d.f. is normal, (2) involves powers of the cumulative normal 
and this is unpleasant to work with. One can avoid this for large n, using a 
technique given in Cramér [2; p. 371]. 

Let us define the random variable 7, as 


(12) tm = nF(y,). 
Then it follows, for any fixed u in 0 < u < n, that 


(12) Prob ( < w) = Prob (nF(y,) <u) = Prob (y, < F“(Y)). 


Substituting in (2), it follows easily that the distribution function of », is given 
by 


(14) Tu) =1-— (1 a “)" 


As n — o, the sequence of random variables 7, converges in distribution 
to a random variable 7 since tne sequence of distribution functions T',(w) con- 
verges for any u to a distribution function, 

15) T(u) = lim T,(@u) = 1 —e™, u = 0. 


n7© 


From (15) it follows that the density function of 7 is given by 
(16) y(u) = lim y,(u) = e™, u > 0. 


no 


It is clear from (12) that the sequence of random variables Ya converges in 
distribution to a random variable y where 


a7 v- 


Hence if we can compute F~‘(n) either exactly or approximately, then we can 
find the distribution of y. 
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Example 1: Suppose that the underlying distribution of X is uniform in (0, A); 
then (12) becomes 


t™m = ny,/A. 
This implies that 


<= 
Yn n n- 


In other words, y, = min {x, , 22, +++ , a} is for large n distributed exponentially 
with the cumulative distribution function, 


r@=1-+6™, v2 
= 0, elsewhere 


and density function 
N ny 
wy) =Zew*, y20 


= 0, elsewhere. 


Example 2: Suppose that the underlying distribution of X is exponential; 
then (12) becomes 


dun 


nl —e ) = ™- 


Hence 


1-- 
nN, 


wie lo : 
Yn Xx g n 


Consequently, to within terms which are o(2), 


Yn ~ 


mr 
Thus y, = min {2 , 22, --- , Za} is for large n distributed exponentially with 
the cumulative distribution function, 
rLQ=-1-«™, y20 
0, elsewhere 
and density function 
yy) =nre™, y>O 
= 0, elsewhere. 


Actually we know, from Equations (10) and (11), that y, is exactly distributed 
according to I,(y) and y,(y) for all n. 


Example 3: Suppose that the underlying distribution of X is described by 
the probability density function 
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f{@) =ar*", O<2z<1, a>O 
= 0, elsewhere. 
In this case, 


F(a) 


It follows from (12) that 


Hence 


From this it follows that for large n, y, = min {2z, , %2, «++ , 2} is distributed 
according to the distribution function 


rg) =1-e"", y20 
= 0, elsewhere 


and density function 


a-1 —ny@ 


Ya(Y) = nay é ’ y=>0 


= 0, elsewhere. 


Remark: Gumbel [7] calls limiting distributions of the kind obtained in Ex- 
amples 1, 2, and 3, type 3 asymptotic distributions of smallest values. The 
distributions obtained in Example 3 occur in statistical theories of strength 
[8] and are frequently called Weibull distributions. Note further that plotting 
log log [1/(1 — I.,.(y))] against log y gives a straight line. Type 3 asymptotic 
distributions arise when two conditions are met: (i) the range over which the 
underlying density function is defined is bounded from below, i.e., F(z) = 0, 
x < 2 for some finite xz , and (ii) F(x) behaves like B(x — x,)*, for some a, 
B > Oasxz— 24. In that case, y, = min {z, , 2, --- , 2} is distributed as- 
ymptotically like the random variable x) + (n/n8)”*. The associated cumula- 
tive distribution is 


r@-1-r~", sem 


= 0, yY< 2%. 


Example 4: The underlying distribution of X is described by the density 
function 
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It follows from (12) that 


Hence 


From this it follows that for large n, y, = min {z, , 2, °-: , X,} is distributed 
according to the distribution function 


r.y =1—e”", y<0 
ol, 29 
and density function 


non 
¥n(Y) - ye y< 0 


= 0, y => 0. 


Remark: Gumbel [7] calls a limiting distribution of the kind obtained in 
Example 4 a type 2 asymptotic distribution of smallest values. This type of 
asymptotic distribution (also called Cauchy) arises when the range of the 
underlying density function is unlimited from below and if for some positive 
k and A, lim,.... (—2)"F(z) = A. In that case the random variable 

Yn = min {x >to, #43 » Tn} 
is distributed asymptotically like the random variable —(nA/n)'*. It can be 
shown that the associated cumulative distribution function is given by: 
ry) =1—-—e"—™" oy <0 


i. y= 0. 


Example 5: The underlying distribution of X is described by the density 
function 


f@)=43e"', -2 <4<e., 
In this case, 
F(z) = }e’, -o<2<0 
1—%4e*, 2x>0. 
It follows from (12) that 
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(this is reasonable since for large n it is virtually certain that 


yo = min {2, ,%,°°* , Ze} 


will lie in the region —~ < y < 0). Thus we get 


Yn ~ —log (3) + log n. 


From this it follows that y, = min {z, , 22, +--+ , Z,} is distributed asymptotically 
with the distribution function 


r,.(y) _ 1 is g ame 
and density function 


u~(n/2)e¥ 


n 
Ay) = 3° 


Remark: Gumbel [7] calls a limiting distribution of the kind obtained in 
Example 5 a type 1 asymptotic distribution of smallest values. Note that plotting 
log log [1/(1 — TI..(y))] against y gives a straight line. 


Example 6: The underlying distribution of X is described by the density 
function 


f(z) = een sie ae ete ae, 
T 


In this case 


F(a) = [ Feet 


It follows from (12) that 


Qn = nF (yn) = £ 5 Feet 


In [2], pp. 374-375, it is shown that asymptotically 


——— , log logn + 4x log 9 
2 = —V2 logn =a, + B, log n, 
. —— 2+/2 logn V/2 log n es 


where 


—— , log logn + 4x 1 
= —V2logn + SSR 0 alae 
. i 2/2 logn V2 logn 


From this it follows that y, = min {x , 22, +++ , 2} is distributed asymptotically 
with the distribution function 


r,(y) ~< en P(e and/Bnd 


Remark: This is another instance in which a type 1 asymptotic distribution 
of smallest values is obtained. As in Example 5, log log [1/(1 — I.,(y))] plots 
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as a straight line against y. A type 1 asymptotic distribution occurs if f(z) 
tends to zero exponentially as z — —, e.g., as e”'*!”, In this case the distri- 
bution of y, = min {2 , %, , «++ , %»} is for large n distributed essentially as 
a, + 8, log n. The associated cumulative distribution function is 


r.(y =1- en tzPl (wan) /Bnd 


In Example 5, a, = —log (n/2) and 8, = 1. In Example 6, 


log log n + 4r 
» = —V2 logn 
m oF eee 


1 
= V2 log n 
3. DIsTRIBUTION OF LARGEST VALUES. 


Up to this point we have dealt exclusively with the distribution of y, , the 
smallest value in samples of size n. Let us now briefly give similar asymptotic 
results for z, , the largest value in samples of size n. The distribution of z, is 
obtained as follows: 


(18) H,(2) = Prob (, < z) = Prob (all z; < 2) 
= Prob (x, < z, 42 <2, °** »%m <2) 
= [F@!. 


H(z) is the cumulative distribution function of z, . The density function h,(z) 
is found by differentiating H,(z). This gives us. 


(19) ha(2) = Hi@) = nf@(F@T. 


Let us consider special cases for which we can get an exact distribution of z, 
for any n. 


Example 1: Underlying distribution is uniform over [0, A]. In this case 
F(z) =0, «<0 
x 
a O<2+<A 
=1 2> A. 


H,{z) = 0, 





ELEMENTS OF THE THEORY OF EXTREME VALUES 


nz” 


A? 
= 0, elsewhere. 
Example 2: Underlying distribution is exponential. In this case 
F(x) = 0, z<0 
=l1-e™, 2z#>0. 
H,@) = [F@r=(l-e™)", 220 


= 0, elsewhere 


h,(2) = Hi) = 0<z<A 


h,(2) = Hiz) =nre “(1 —e™)*", 22>0 
= 0, elsewhere. 


In general it is not easy to write down an explicit and (or) tractable formula 
for H,(z) or h,(z). The following method [2; p. 371] is very useful for getting 
asymptotic formulae for large n. 

Let us define the random variable é, as 


To find A,(v), the distribution function of  , we note that for any vin0 < v <n, 


(21) A,(v) = Prob (, < v) = Prob (Fe) >1- 2) 


Prob (. > r( ~ 2) =1- nf (1 _ 2) ] 
[He 2) 1 G3 


As n — o, A,(v) converges for any v to a distribution function A(v) where 
(22) Av) = limA,@) =1-—-e", v20. 


Thus the random variables ~, converge in distribution to a random variable é. 
The density function of ¢ is given by 


(23) \v) = lim>,@) =e", v>O0. 


It is clear from (20) that the sequence of random variables z, converges in 
distribution to the random variable z, where 


a -r(i-2) 
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Hence if we can compute F-‘(1 — £/n) either exactly or approximately, then 
we can find the distribution of z. 


Example 1: Suppose that the underlying distribution of X is uniform in 
(0, A); then (20) becomes 
or) 


“Z,~ A 4. 
n 


This means that the distribution of largest values can, in this case, be expressed 
as the difference between a constant (the upper limit of the range of possible 
values) and numbers drawn at random from an exponential distribution. 
Analytically the cumulative distribution function can be written as 


= 0, z2<0 
A,(2) a on 0 < z < A 
1, z>A 


d,(2) = Al(2) = ier 0<z2<A 


= 0, elsewhere. 


Remark: Gumbel calls limiting distributions of the kind obtained in Example 
1, type 3 asymptotic distributions of largest values. Type 3 asymptotic distri- 
butions of largest values arise when two conditions are met: (i) the range over 
which the underlying distribution is defined is bounded from above, that is, 
F(x) = 1, x > po, for some finite z and (ii) 1 — F(x) behaves like B(2. — x)*, 
for some a, 8 > 0 as x — a . Under these conditions z, = max {x, , 12, --* , Zn} 
is distributed asymptotically like the random variables 1 — (&/n@)*. The 
associated cumulative distribution function is 


A,(z) a eee zZ < Lo 
= 1, 2 Ei 


Example 2: Suppose that the underlying distribution of X is described by 
the following density function 


f(z) 
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In this case 


&, = n(1 — F@,)) = = 


n 


eo a. Se 


g 


From this it follows that for large n, z, = max {x, , 42, --+ , 2} is distributed 
according to the distribution function 


A,(z) = 0, z<1l, 
a . : > 1 


and density function 


A(z) = ee z2>1 


= 0, elsewhere. 


Example 3: Suppose that the underlying distribution of X is described by 
the following density function 


f(x) = 0, 


In this case 


& = n(l — F@,)) = n/en 


(2) 
oc Zn ™ ¢ ° 


From this it follows that for large n, z, = max {z, , 22, -++ , Za} is distributed 
according to the distribution function 


A,(2) = 


and density function 


= 0, elsewhere. 
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Remark: Gumbel calls limiting distributions of the kind obtained in Examples 
2 and 3, type 2 asymptotic distributions of largest values. This type of asymp- 
totic distribution occurs when the range of the underlying density function is 
unlimited from above and if for some positive k and A, lim,.. z*(1 — F(z)) = A. 
It can then be shown that z, = max {z, , 72, -*+ , a} is for large n distributed 
asymptotically like the random variable (nA/t)”*. From this it follows that 
the associated cumulative distribution function is given by 







A(z) = 0, for 2<0 














—nA/2* 
? 


=e for 22> 0. 
In Example 2, k = 1 and A = 1; in Example 3, k = aandA = 1. 


Example 4: The underlying distribution of X is exponential with 







F@=0, «<0 


In this case 


& = n(1 — F(,)) = ne 


and therefore 







2 ~ em _ loge, 
7 A A 










From this it follows that for large n, z, = max {2 , 22, - 


- , 2} is distributed 
according to the distribution function 








A,(2) _— 0, ‘<= 0 


——". ens 










and density function 
rm) tt A‘) = ner 


= 0, elsewhere. 


Remark: Note that in this case a double exponential density of the form 


pager 


é 


(we neglect constants) is involved. Gumbel calls such distributions type 1 
asymptotic distributions of largest values. Such asymptotic distributions 
are obtained if the underlying density function falls off exponentially as s — © 
(e.g., f() ~ e74* as x — &). When this happens, z, = max {% , Ta, °°: , Ze} 
is distributed for large n as a, — 8, log & (where 8, > 0). For example, if 
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—-o <Z< @, 


7 ~~ _ log logn + log 4x 
~= 721 -— OO 
oo = aes 2+/2 logn 


1 
9/2 log n 
In the pure exponential case just treated 


Bn 


— logn 
a, = x 
and 8, = 1/d. In distributions of this kind log log 1/A(z) should plot linearly 
against z. 

We have shown how one can find the c.d.f. and p.d.f. of smallest or largest 
values in samples of size » drawn from a distribution described by the p.d.f. 
f(x) and c.d.f. F(x). It is particularly interesting to. find the mode of g,(y) (or 
the most probable value of the smallest value in samples of size n) and the 
mode of h,(z) (or the most probable value of the largest value in samples of 
size n). We shall limit ourselves to the problem of finding the mode of h,(z), 
leaving it to the reader to carry out analogous considerations for g,(y). 

If f(x) is differentiable, then the modal value z* of h,(z) is the solution to 
hi(z) = 0. Thus in many cases this amounts to finding z* such that 


(25) f'@)F@*) + @ — Ife) = 0. 
Thus one will find in Example 2 that 


n+1 
2 ? 


a= 


in Example 3 that 


2? = (mi ta)" 
7 a+l : 


and in Example 4 that 
2+ — loen, 
7 A 


In Example 1 one cannot obtain the modal value by means of (25), due to 
the fact that the density is zero for z > A. In this case 2* = A,.i.e., the right- 
hand endpoint of the interval is the most likely value for the largest value in 
a sample of size n. If the underlying density is normal, there are difficulties 
in solving (25) explicitly. For n large the modal value is given by 


log log n + log 4x 
z=a,= V2) a eR 
" — 2+/2 logn 
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For large n one can also find the modal value by finding the maximum of \,(z). 
The maxima of X,(z) are as follows: 


Ex. 1 


Ex. 2 
Ex. 3 


Ex. 4 


_ log log n + log 4x 
2+/2 log n 


Remark: One can learn a great deal from these examples. Essentially one 
can assert the following: 
(i) If f(z) = 0, for x > A, then the modal value of the distribution of largest 
values will be A. 
(ii) If 1 — F(x) — 0 in such a way that lim,... x*(1 — F(x)) = const., then 
the modal value of \,(z) can be expressed for large n as a + bn”°. 
(iii) If f(x) ~ e~4* as x — @, then the modal value can be expressed for large 

1/k 


nasa + b (log n)””. 


Normal distribution 4/2 logn 


4. Remarks On THE AppLicaTIONs Or EXTREME VALUE THEORY. 


The primary aim of the paper is to present some of the basic elements of 
extreme value theory. It should be pointed out that the theory can be applied 
in many ways. The relevance of this theory to floods, droughts, and other 
meteorological phenomena, such as extreme pressures, temperatures, etc., is 
clear. Still other examples are furnished by fracturing or fatigue of metals 
and textiles, by the breakdown of dielectrics and the corrosion of metals. In 
the case of fracturing or fatigue what may be involved is essentially a break- 
down of the weakest link in a chain; in the case of dielectrics, the breakdown 
will generally occur at the spot where the dielectric has been penetrated most 
deeply by defects known as conducting particles; in the case of corrosion, pene- 
tration of the metal may well be associated with what happens at the largest 
pit. 

What does the statistical theory of extreme values contribute toward the 
solution of the physical problems that we have just discussed? Two important 
contributions are: (i) It explains in a theoretical way why certain distributions 
may be expected to occur, and this makes it easier to analyze extreme value 
data analytically or graphically. For example the consideration in Section 
3 helps to explain why one plots log log 1/A(z) versus z in dealing with meteoro- 
logical extremes. (ii) One can obtain a rational explanation of the size effect 
(this is the name given to the phenomenon that larger specimens will, for example, 
generally fracture with less applied stress, break down with less applied voltage, 





ELEMENTS OF THE THEORY OF EXTREME VALUES 41 


or corrode in a shorter time). For if the behavior of a specimen depends on 
the smallest or largest value (as the case may be) of some defect or flaw, then 
increasing the specimen size may mean that we are taking extremes from larger 
and larger samples. Hence the modal value should change. We have seen in 
Section 3 that under certain conditions the modal value of the largest [smallest] 
value in samples of size n changes as a + bn”“[a — bn'’“], under other con- 
ditions as a + b (log n)““[a — b (log n)*], and in still other cases is independent 
of n. This can be immediately interpreted as a dependence on volume or area 
or length with n replaced by V (volume), or A (area), or L (length). 

More specifically, it is shown in [4] and [5] that the theory of extreme values 
is relevant to the description of the distribution of breaking strengths and 
of the size effect. In [6] extreme value theory is applied to both the distribution 
and size effect problem in connection with the dielectric strength of paper 
capacitors. In [1] it is shown that maximum pit depth data taken from aluminum 
samples follow an extreme value distribution of the doubly exponential form 
and in [3] it is reported that maximum pit depths in areas of size A increase 
linearly with log A. 

Frequently the statistical study of engineering data is greatly facilitated by 
plotting the data on appropriate graph paper. Gumbel has constructed special 
graph paper for the plotting of extreme value data. Many examples of the 
use of such paper in analyzing meteorological and fatigue data are given in [7]. 
Graphical analyses of corrosion data are given in [1] and [3]. 
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System Efficiency and Reliability* 


R. E. Bartow Anp L. C. HUNTER 


Sylvania Electronic Defense Laboratory 
Mountain View, California 


This paper presents a method for determining the reliability of large, complex 
systems. Repair is an integral part of the model proposed and the usual assumption 
of component independence need not be made. Recognizing the stochastic nature 
of an electronic system operating in time, the reliability of any system is defined 
in terms of a suitable stochastic process. This enables the use of powerful proba- 
bilistic techniques which have been developed in recent years. In particular, the 
theory of Markov processes can be readily applied. Given component repair and 
failure parameters, the reliability of the system at any specified time as well as the 
limiting reliability can easily be obtained. If the reliability requirement for the 
system is given, it is a routine matter to determine permissible values for the mean 
time to repair and the mean time to failure of each component. 


INTRODUCTION 


Most analytical papers in the field of reliability have three major drawbacks: 
(a) Components are assumed to operate independently of one another. (b) Only 
the mean time between failures (MTBF) is investigated. (c) Repair is not an 
integral part of the reliability model. In reality, however, the electronic system 
can be repaired, the components are not independent, and one is interested 
not in the length of time between failures, but in the probability the equipment 
will be operating at any given future time. 

In this paper we are interested in the general problem of assigning a reliability 
function to a system which overcomes these three drawbacks. To do this we 
must abstract the system idea to a mathematical model. To succeed, this model 
must be able to cope with the following handicaps: (a) Components of large 
systems rarely operate independently. (b) Meager information is usually avail- 
able concerning failure and repair parameters. The model must be sufficiently 
complex to allow for the first handicap and sufficiently simple to cope with the 
second. In particular, the model should be sufficiently simple to allow a study 
of the effects of varying the input parameters over a wide range. 

A brief background concerning the statistical techniques used in this paper 
is perhaps appropriate. Many authors have considered the effect of repair 
on systems consisting largely of identical, independent operating units. His- 
torically, C. Palm [6], following the methods of A. K. Erlang [1], was one of 
the first to use the properties of the exponential distribution in solving various 


* This paper was presented at the IRE National Convention, 1959, New York City and 
was prepared under Signal Corps Contract DA-36039 SC-78281. 
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“repairman” type problems. Feller [3], popularized and extended the technique. 
More recently, renewal theory techniques have been used by many authors 
in solving related problems in telephone trunking, inventory theory, and elec- 
tronic counters [7]. A good discussion of definitions for reliability work is given 
by Knight, Jervis and Herd [4]. However, no natural extension of the definition 
of reliability to include repair appears. 

In the next section we extend the usual definition of system _—— and 
reliability to include repair. This is followed by an application to a specific 
system and a summary. 


MatTuHEMATICAL MopELs For System 
EFFICIENCY AND RELIABILITY 


First of all, we note that a system can be in one of possibly many states at time 
t. To take a most elementary example, the states could be the ‘‘on” state (Ep) 
and “‘off’’ state (Z,). In general, the state space will be determined by the physical 
configuration of the system as well as the raison d’etre for the system. We use 
S to denote the (possibly uncountable) state space of the system and E; to denote 
a particular state. 

We require some method for describing the manner in which the states are 
realized in time. It is apparent that the state of the system at any time ¢ is really 
a random variable. For convenience we will take 0 to be the origin of the time 
parameter. These considerations lead us to the following definition. 

Definition 1. The state of a system as a function of time is a stochastic process, 
which we denote by {X(#), 0< t < ©}, measurable with respect to the usual 
background probability space (Q, A, P). 


A sample function of the process, that is to say any particular realization of 
the process, will be denoted by z(/). As an example, consider the following 
realization of a one-unit system capable of assuming only two states (on-off). 
The step function is our x(#). Note that the 


t, Ms 
Fig. 1—Realization of a Two State System. 
system was originally on, that it failed at ¢, , was repaired at time /, , etc. 


As a simple example of a system with a continuous state space let S be the 
space of all possible voltage levels corresponding to the voltage output of a 
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linear vacuum tube circuit. Due to a phenomenon called the “shot effect”? we 
obtain a fluctuation of the voltage level about a mean value in time. 

In general, the state space describes possible outcomes for the system while 
the stochastic process tells how the system can operate in time. The system 
changes from one state to another when a component fails or when a component 
is repaired. In practice each component is usually considered to be in one of 
two possible states (on-off). In this case if there are n components, the number 
of possible states for the system is N = 2”. Fortunately, other considerations 
usually reduced the number of states considerably. 

Usually it is more advantageous for the system to be in a certain class of 
states which we may call “favorable”. The designation of these favorable states 
and the choice of appropriate weights for the states in the state space leads us 
to define another element of our model. For example, consider again the basic 
two-state system. If the system is on at time ¢ we could consider this a unit 
gain or payoff. If the system is off at time ¢ we could consider this to be either no 
gain or no payoff. This leads us to define a gain function or payoff function. The 
gain function should assume values on the real line so that we can evaluate our 
gain quantitatively. Hence, as a first approximation it is clear that the gain 
function, designated g, is a Borel measurable function defined on the state space 
S of the system. Since the state of the system as a function of time is a stochastic 
process the gain function should be composed with the process, 7.¢., g[X(2)]. 
This now defines a new process superimposed on the old. As a simple example, 
suppose that the system is a single vacuum tube. Furthermore, suppose that its 
time-to-failure is distributed exponentially with mean life 1/A. Let us consider 
it a gain of one unit if the tube is burning at time ¢ and zero otherwise. One 
realization of g[X(t)] might appear as follows: 


gLx(t)] 


‘, 


Fig. 2—Realization of a Gain Function. 


Since we wish to describe the reliability of the system as a single function 
of time, obviously we cannot use g[X(¢)], however the expected value of this 
process can be used. 


Definition 2. The reliability of a system is 


RW) = F(X = f olXl, w)] aPw). 





46 R. E. BARLOW AND L. C. HUNTER 


This can be phrased as the reliability is the expected gain of the system at time /. 
To illustrate this definition of reliability, let A denote a class of favorable 
or acceptable states. Then we might define 


aco) = | if X(teA 
0 otherwise 


Thus, R(¢) would be the probability that the system is in one of the favorable 
states at time ¢. For the single vacuum tube system R(t) = E{g[X()]} = 
exp (—?#/A). The graph of R(#) appears as follows: 


R(t) 


Fig. 3—Reliability of a Single Vacuum Tube. 


One other factor must be considered. A system may be operating against a 
threat which is more likely to appear at certain times than others. The environ- 
ment may provide a distribution function F on the time axis. For example, 
the threat might appear at any time in the interval [0, 7] with equal probability. 
Then 


if ¢<0 
F,() = if O<t<T 


if otherwise 


This leads us to attach a numerical value to the system independent of the 
time parameter which we will call the efficiency of the system. 


Definition 3. The efficiency of the system is 


m= [ BtolX(On ar 
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Note that if F = F, then Eff is merely a time average of the reliability of the 
system in the interval [0, 7]. In ge=cral, the efficiency can be phrased as “‘the 
expected reliability of the system.” If we confine ourselves to time averages then 


Effe = pf BlolXN} at 


It will be convenient to define the limiting efficiency of a system in terms of 
a time average. 


Definition 4. The limiting efficiency of the system is 
Eff. = lim Effr 


For the single vacuum tube system, using the g defined before and letting 
F = F, we obtain 


Effe = pf exp (—t/0) dt = dU — exp (-T/I/T. 


As T approaches infinity the efficiency approaches zero since “most” of the 
time the vacuum tube will be dead. Note also that R() = 0. For this example 
Eff, is the expected fractional amount of time that the tube is burning in the 
time interval [0, 7’. 


A more interesting example is the one-unit two-state system described earlier. 
Let 


1 if X(é) = FE, 


g{X(t)] = 
0 if X() =£,. 


Suppose that the system fails according to an exponential distribution with 
mean 1/) and is repaired according to an exponential distribution with mean 


1/u. Then the reliability of the system is the probability that it is on at time ¢ 
which is 


— r mm 
It is easy to show that the expected fractional amount of time that the system 
is on in the interval (0, 7] is 


a r _ AXexp[-—Q + wT] 
A+u TA+p)” TA + ») 


The above definitions have been made using an underlying probability space 
(2, A, P). However, in order that these definitions be readily applicable, it 
should be possible to evaluate system reliability, efficiency, etc., without knowing 
either P or &. In practice, we only know the various components of the system 
and something about the manner in which they fail. Fortunately, it is usually 
unnecessary to determine 2 and P. 

The model presented is very general, indeed, and it must be admitted that 


Effr o 
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actual computation of system efficiency and reliability is easily accomplished 
only for systems with finite state spaces and for those which can be described 
by Markov processes. Recall a Markov process is a stochastic process such that, 
given the system is in state EZ; at time & , no additional information concerning 
states of the system at previous times can alter the (conditional) probability 
of the process being in state E; at some future time ¢. Expressed in another way, 
the future states of the process depend only on the immediate past history 
of the system. The reader is referred to Doob [2], and Feller [3], for a complete 
discussion of Markov processes. 

A few qualitative statements should perhaps be made before an example 
is given. Any system whose individual components fail approximately according 
to an exponential distribution is best described by a stationary Markov process. 
If we want repair as a feature of our model, stationarity is essential to avoid 
extremely complicated analyses. The system is best described by a non-stationary 
process if the failure rates of components change with time. 


Examp_e: A THREE Unit System 


Consider a simple system consisting of two transmitters and power supply 
providing a high voltage with modulator and modulation programmer. We are 
interested in the probability that this system will be operating properly at any 
specified future time, which for this example will be the “reliability” of the 
system. Since even this simple system could consist of many non-independent 
components, the answer to this question is not apparent. 

We begin by considering the following three units which compose the system: 
Transmitter A, transmitter B, and a main unit consisting of the power supply. 
Figure 4 presents a system block diagram and the estimated failure and repair 


TRANSMITTER A TRANSMITTER 8B 


\/A = 75 HRS IF TRANSMITTER B A = 100 HRS IF TRANSMITTER A 
IS OPERATING IS OPERATING 

4 = 50 HRS IF TRANSMITTER B A = 50 HRS IF TRANSMITTER A 
IS NOT OPERATING IS NOT OPERATING 

Wp = 2 HRS ih = 2 HRS 


= MEAN TIME TO FAILURE 
= MEAN TIME FOR REPAIR 
OR 
FAILURE RATE PER UNIT TIME 
= REPAIR RATE PER UNIT TIME 


MAIN UNIT 


A = 200 HRS 
4 = 5 HRS 


. 4—A Three Unit System. 
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times of each unit. Note that the mean time to failure of transmitter A is 75 
hours if transmitter B is operating; 50 hours if transmitter B is not operating. 
The mean time to failure of transmitter B is 100 hours if transmitter A is oper- 
ating; 50 hours if A is not operating. This hypothetical system is representative 
of transmitters using mechanically tuned magnetrons. Two transmitters are 
used, each tuning one half the desired frequency range. However, if one trans- 
mitter should fail, the other can tune the entire range with a resultant decrease 
in the expected time to failure. 

At any specified time all three units could be operating properly. In this 
case we shall say the system is in state EZ, . However, there are other possible 
states for the system. We shall say the system is in state 


E, if transmitter B and the main unit are functioning properly and trans- 
mitter A is inoperative, 

E, if both transmitters are working properly and the main unit has failed, 

E, if transmitter A and the main unit are working properly and transmitter 
B has failed, 


E, if both transmitter A and the main unit have failed and transmitter B 
is functioning correctly, 
E, if both transmitters have failed but the main unit still functions correctly, 


E, if transmitter A, is working but transmitter B and the main unit have 
failed. 


If the main unit fails the system as a whole becomes useless and must be 
shut down while repaired. Similarly, if both transmitters fail, the system serves 
no purpose and operation is discontinued until repaired. Hence, the other possible 
state for the system, where all units have failed, is ignored. 


B AND MAIN A ANDB A AND MAIN 
UNIT FAIL FAIL UNIT FAIL 


E € 
7 oe , 
has As Mg 


° 
h36 


y v ul 
~~ 35 MAIN _ FAILS 'S ‘a 
3 2 i 

rae 


02 Ao, A FAILS 


B FAILS 


Hos Hoo Ho 
v 
™! ra 
. 5—State Space for a Three Unit System. 
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System reliability and efficiency at time ¢ can be specified if we know the 
probabiity we are in each of the states E, , E, , --- , E, at time ¢. To find these 
probabilities, we assume the time-to-failure measured from the time of the 
last repair, and the time-for-repair of each component are exponentially distri- 
buted random variables. The Greek letter \ with subscripts will always be used 
to denote failure rates per unit time. The Greek letter » with subscripts will 
denote repair rates per unit time. Figure 5 represents the state space for the 
system with the associated distribution parameters. For example, Ao; is the 
distribution parameter in going from state E, to state E, . In other words, 
1/Xo, is the mean time-to-failure of transmitter A when all units are operating 
correctly. 

Arrows in Figure 5 denote directions to possible neighboring states. We have 
assumed for this example that whenever the main unit fails, all available repair 
effort is devoted to it. Hence, transitions from E, to E, and E, to E, are impossible. 

Let P,(#) denote the probability of being in state EH, at time ¢ given that 
rr 0(0) = 1, 

Referring to Figure 5 for the distribution parameters, we can write the fol- 
lowing system of equations: 


Po(t + h) = Po(t)[1 — dosh — Yosh — Nosh] + P,(A)uosh 
+ Pa()uosh + Pa(duosh + oh) 
P,(t + h) = P,(t)[1 — woh — Arh — Arsh] + Po(t)rork 
+ P,(i)ussh + Ps(t)uish + o(h) 
P(t + h) = P2(t)[1 — mock] + Polt)Aceh + o(h) 
P,(t + h) = P,(1)[1 — posh — Asch — rash] + Po(t)roah 
+ Pe(t)usch + Ps(t)ussh + o(h) 
Pit + h) = Pal — wash — wish] + Pi()Auh + o(h) 
P,(t + h) = Py(f)[1 — mash — mish] + Pa(t)rash + Pi(é)Arsh + o(h) 
P,(t + h) = P,(t)[1 — sch — posh] + Pa(t)rsch + 0(h) 


We illustrate the line of reasoning used to arrive at the above equations by 
deriving P,(¢ + h). If the system is in state E, at time ¢ + h, then five possi- 
bilities arise: 


(a) The system was in state EZ, at time ¢ and no change occurred. 

(b) The system was in state EZ, at time ¢ and a transition to Ey occurred 
during (f, ¢ + h). 

(c) The system was in state EZ, at time ¢ and a transition to EZ, occurred 
during (¢, ¢ + h). 

(d) The system was in state FE, at time ¢ and a transition to Ey occurred 
during (¢, ¢ + h). 

(e) The system is in state E, , E, , or E, at time ¢ and state E, at time ¢ + h. 
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According to our hypothesis of exponential failure and repair time, the proba- 
bility of the first contingency is P,(#) times the probability of no change during 
(t, ¢ + h) and this last is 1 — Aoih — Aosh — Aosh + o(h). Similarly, the second 
contingency has probability P.(t)uooh + o(h) while the third and fourth coa- 
tingencies have probability P.()uo2h + o(h) and P3(t)uosh + o(h), respectively. 
The fifth possibility has probability of smaller order of magnitude than A. 
Summing these probabilities gives Po(¢ + h). 


Letting h — 0 in each of the equations, we obtain the following system of 
linear differential equations: 


Po(t) = —(Aor + Ao2 + Aos)Polt) + Pi()uor + Po(t)uoe + Pa(é)uos 
Pi(t) = —(uor + Ara + Ars)Pi(f) H+ AorPolt) + wraPa() + masPs(d) 
P(t) = —Ho2Po(t) + do2Po(t) + maePa(t) + peePol(t) 

P3(t) = —(uos + Ass + Ase)Pa(t) + AosPolt) + wsePo(t) + ussPs(2) 
Pit) = —(Hoa + wra)Pa(t) + AuPi(O) 

P(t) = —(uss + wis)Ps(t) + AssPa(t) + ArsP1 (2) 

Pot) = —(use + uae)Po(t) + AsoPa(d) 


(Note: The above system of differential equations is merely the Chapman 
Kolmagorov equations for stationary Markov processes and could have been . 
written immediately. However, in practice, it is often easier to employ an argu- 
ment of the type presented. See [3].) 

Referring to Figure 4 for the failure and repair rates per unit time, 

Pi(é) = .0133P,(é) — .525P,(#) + .2P,(é) + .5P,(d) 
Pit) = .005P,(é) — .2P.(2) 

P(t) = .O1P,(t) — .525P3(é) + .5P;(t) + .2P,(é) 
Pi(t) = .005P,(t) — .2P,(d) 

P(t) a .02P,(%) + .02P;(é) a P,{t) 


These equations are easily solved by standard methods. Perhaps the easiest 
and fastest method of all is to compute the solution using an analogue computer. 
Table 1 gives the values of P,(¢), P,(¢), P2(é), and P3(¢) when ¢ = 0, 1, 2, 4, 8, 16. 


Time in Hours 


t=1 t=2 t=4 t=8 t = 16 
.9774 - 9630 9475 - 9367 - 9323 
.0102 -0162 -0217 0243 .0247 
.0045 -0081 .0133 -0189 .0224 
.0077 .0122 -0164 .0185 -0187 
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The reliability of the system at time ¢ will in general be a weighted average 
of P,(¢), P,(é) and P;(¢). Using equal weights we may define the reliability of 
the system as 


Ri) = Pot) + Pi) + P3(d) 


The limiting reliability of the system is particularly easy to calculate. Recall 
ast— o, P-(t) +0 and P,(t) — P, . Hence, letting t + © equations (2) become: 


— .0283P, + .5P,.+ .2P.+ .5P; =0 
-9133P, — .525P, + .2P, + .5P; = 0 
.005P, — .2P, =0 
01P,) — .525P; + .56P; + .2P, =0 
.005P, — .2P, _ =0 
.02P, + .02P, — P; =0 
.005P; — .02P, =0 

The solutions to this set of algebraic equations are: 
Py = .9314 
P, = .0247 
P, = .0233 
P; = .0187 
P, = .0006 
Ps = 0009 
P, = .0004 


Hence the limiting reliability of our system is 


R(e) = Pp +P, + P; = .9748. 


SUMMARY 


It has been noted that an electronic system can be in one of many possible 
states and that state of the system at any specified time is a random variable. 
This leads us to consider the state of a system as a function of time to be a 
stochastic process, and to attach a value to each possible state. This value 
measures the payoff or gain to the system’s user when the system is in that 
state. The reliability of a system at any specified time was defined as the expected 
payoff or gain at that specified time. The efficiency of a system over a time 
interval [0, 7] was defined as the time average of the reliability over this time 
interval. It was noted that the usual definitions of reliability and efficiency 
are special cases of these more general definitions. 

Given a system for which a reliability analysis is to be made, the first task is 
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to break the system down into its primary components and determine the 
state space. An example was given of how this can be accomplished. In parti- 
cular, if the waiting times between changes in the states of the system is assumed 
to be exponential, the operation of the system can be described by a Markov 
process and the familiar techniques of Feller [3] employed. 
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This paper presents a new table and some approximating polynomials especially 
designed to facilitate maximum likelihood estimation of the parameters of the 
gamma distribution, and also applicable to the 2-parameter Type V; it discusses 
methods of computing the sampling variances of the likelihood estimators; it illus- 
trates the use of the tables in a numerical example; it mentions applications to the 
Erlang distribution and difficulties of application to the general Type III. Finally 
it inquires when the numbers extracted from the tables are maximum likelihood 
estimates, and what they are estimates of. 


I. TuHeory: LIke.tinoop Estmators For Tue Type III 


“Curves of Pearson’s Type III. offer a good example for the calculation of 
the efficiency of the Method of Moments.” R. A. Fisher (1922, p. 332) then 
went on to derive expressions for maximum likelihood estimators of the parame- 
ters of the Type III, and thence to compute the variances and efficiencies of 
the corresponding moment estimators. He concluded that the method of moments 
is inefficient, except for distributions closely resembling the normal. 

Although curve fitting has lost its general popularity, the Type III, and 
in particular the 2-parameter gamma distribu‘ion, continues to be fitted. The 
general Type III density function may be written: 


1 — m\*" 

al(p) (? a ) exp [-—(@ a m)/a). (1) 
According to Fisher’s derivation (cf. Kendall, 1946, example 17.19), the logarithm 
of the likelihood function for a sample of n observations is (with some rear- 
rangement and with p = p — 1 in Fisher’s notation) 


= —nplna —nlInI(p) + (p — 1) DS In(@e — m) — D(x — m)/a. (2) 


* This is the revision of a paper read before the meeting of the Institute of Mathematical 
Statistics, Cambridge, Massachusetts, August 28, 1958. 

The underlying research was supported in part by the U. S. Army Ballistic Missile Agency 
under Contract No. DA-11-022-ORD-2732, and in part by the Sloan Research Fund of the 
School of Industrial Management at the Massachusetts Institute of Technology. 

t Now at Princeton University. 
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Then, the maximum likelihood estimators 4, 6, and m are the roots of the three 
simultaneous equations 


nap = >. (x — m) (3a) 


—n +-In To) + Dn @ — m) = nina (3b) 


1 
== 


a(p — 1) 


obtained by differentiating (2). The solution for all three unknowns is cumber- 
some: moreover, it encounters grave difficulties when p < 1, i.e., for J-shaped 
curves; equation (3c) cannot then be satisfied, and m must be estimated from 
the smallest observation. 

For the gamma distribution, when m is known, the simultaneous solution of 
(3a) and (3b) for a and p is fairly convenient and always possible. If A represents 
the arithmetic mean of the deviations x — m, and if G represents their geometric 
mean, the substitution of a = A/p from (3a) into (3b) gives 


oe (30) 


n=y (4) 
where 


d 
quinp—7 mri): 


y=nA—1n4%, 


Solution of (4) with the aid of a table of natural logarithms and a table of the 
digamma function is too laborious for routine use. 


II. TABULATED SoLutions AND THEIR UsE 


Including those presented herewith, the following tables are available to aid 
in the solution of (4): 
Function Description Reference 


n(p) 5D p = 0.01(0.01)2(0.02)5(0.1) Chapman* 
20(1)100 


n(p) 7D p = 0.1(0.05)3(0.1)5(0.5)10 Masuyama and 
(5)20(10)50 Kuroiwa (1951) 
n(p) logice 7D same game 
e” 8S same 
pn(p) logice 7D same same 


p = ¢(n) 2D 7» = 0.01(0.001) 0.1(0.005) Chapman (1956) 
0.5(0.01)0.57 


38D 7» = 0.58(0.01)1 same 
np = n(n) 8D » = 0.00(0.01)1.4 Table 1A 
7D » = 1.4(0.2)18 Table 1B 


* Chapman (1956) reports that this table is available in mimeographed form from the 
Laboratory of Statistical Research, University of Washington. 
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Tables of n(p) have two disadvantages: they require inverse interpolation, 
and the function varies so rapidly with p that accurate interpolation is difficult. 
Tables of p = $(n) eliminate the first difficulty by permitting direct interpolation. 
Tables of pn(p) eliminate the second difficulty by inducing a function that: varies 


TABLE 1A 
np = plln p — I’(p)/T()] as a function of » = In p — I(p)/T(p). 


” np 6? a 8? 


.00 50000000 —1111 
166108 —1128 | .36 55234633 —1036 | .71 59049158 1.06 62018509 
331088 —1144 | . 358965 —1024 | .72 144048 .07 094026 
494925 —1158 | . 482273 —1013 | .73 238267 .08 169093 
657603 —1171 | . 604568 —1001 |} .74 331825 .09 243714 
819111 —1182 | . 725862 —990) . 424729 -10 317895 


: 979438 —1191 | . 846165 -—978) . 516987 -11 391639 
.07 51138573 —1200 | . 965491 —967) . 608607 12 464951 
.08 296508 —1206 | .43 56083849 -955)| . 699595 .13 537836 
: 453237 —1212 | . 201253 -—944| . 789961 .14 610298 
10 608755 —1216 | . 317712 —933)] . 879710 —609 | 1. 682342 


763057 —1218 | . 433238 —922) . 968850 —602 | 1. 753971 
: 916141 —1220| . 547842 -—910| .82 60057389 —595 | 1. 825190 
.13 52068005 —1220 | . 661536 -—900| . 145333 —588 | 1. 896003 
.14 218649 —1219 | . 774331 —889| . 232689 —581 | 1. 966414 
368074 —1217 | . 886236 -—878| . 319464 —574 | 1.20 63036427 


516282 —1214 | . 997264 -—867| . 405665 —568 | 1. 106046 
663277 —1210 | .52 57107426 —856| . 491298 —561 | 1. 175276 
809061 —1205 | . 216730 —846| . 576369 —555 | 1. 244119 
: 953640 —1200 | . 325189 —836| . 660886 —549 | 1. 312580 
.20 53097019 —1194 | . 432812 —825| . 744854 —542 | 1. 380662 


239205 —1186 | . 539610 -—815| . 828279 —536 | 1. 448370 
380204 —1179 | . 645593 -—805/ . 911168 —530 | 1. 515706 
520025 —1171 | . 750771 —795]| . 993527 —524 | 1. 582675 
658674 —1162 | . 855154 —785| .94 61075362 —519 | 1. 649279 
796161 —1153 | . 958751 —776| .95 156678 —513 | 1. 715523 


; 932495 —1144 | .61 58061573 —766| . 237481 —507 | 1. 781410 
.27 54067685 —1134 | . 163629 —757| . 317777 —502 | 1. 846943 
.28 201741 —1124| . 264927 —748| . 397571 —496 | 1. 912126 
334674 —1113 | . 365479 —738] . 476869 —491 | 1. 976962 
466493 —1103 | . 465291 —729 | 1. 555677 —485 | 1.35 64041454 


597209 —1092 | . 564375 —720} 1. 633999 —480 | 1. 105606 
726833 —1081 | . 662738 —712} 1. 711841 —475 | 1. 169420 
855377 —1070 | . 760389 —703 | 1. 789208 —470 | 1. 232900 
: 982850 —1059 | . 857337 —695 | 1. 866105 —465 | 1. 296049 
.35 55109265 —1047 | . 953591 —686 | 1. 942537 —460 | 1. 358870 
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slowly with p. Finally, our tables of n(n) eliminate both difficulties by permitting 
direct interpolation in a function that varies slowly with 7; hence these tables 
are well suited for the special purpose of solving (4). 

For calculating most of the values in Table 1, we evaluated the function 
np = piln p — d/dp In I(p)] for convenient values of p—using the National 
Bureau of Standards (1953) tables of natural logarithms and H. T. Davis 
(1933) tables of the psi function—i.e., the digamma function; then we obtained 
the inverse function p = ¢(7) by interpolation, which was greatly facilitated 
by the fact that np varies much more slowly than 7 itself. As a convenience 


Tasie 1B 
np = p[ln p — I"(p)/T(e)] as a function of n = In p — T(p)/T(p). 


6? np 3? np 


.8010527 .2 .8536588 
035935 : 548549 
060539 : 560277 
084382 ‘ 571778 
107500 ; 583059 


.6435887 —12935* 
555077 —10841* 
663340 —9204* 
762335 #—7902* 


Noe ee 
CoMOm 
onnas 
CoO fb 


853382 —6852* 
937540 —5993* 
.7015677 —5282* 
088511 —4688* 
156639 —4186* 


129931 594127 
151707 : 604989 
172861 ; 615651 
193420 626118 
213412 636397 


go 89 WON NO 
CoO fw 
90 00 00 oo 
CoO S to 


220568 —3758* 
280727 —3391* 
337487 —3073* 
391165 —2797* 
442040 —2556* 


232862 5 646493 
251795 ; 656411 
270232 : 666156 
288196 j 675733 
305705 . 685146 


3.2 
3.4 
3.6 
3.8 
4.0 


_ 


CM ORN 


490353 —2343* - 322778 694401 
536319 —2155* j 339433 ; 703502 
580126 —1988* ; 355687 ; 712452 
621942 —1839* : 371555 F 721255 
661916 —1705* ; 387052 : 729916 


CoOm & ty 


700182 —1588 ‘ 402193 3 738438 
736861 —1479 : 416989 ; 746825 
772060  —1381 : 431455 ‘ 755080 
805879 —1292 : 445602 ; 763207 
838406 —1211 : 459441 ; 771208 


> or or or or 
CoO » by 


869722  —1137 : 472983 
899901  —1069 ; 486239 
929012 —1007 ; 499219 
957116 —950 ; 511931 
984270 —897 ‘ 524384 


NDABO 
CoO mf te 


* These differences have been modified. 
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and partial check, we obtained a few values by summing special series that 
converge rapidly in the vicinity of zero and one. We carried all computations 
to ten decimal places, although it was clear enough that the tenth place was 
not particularly reliable. Differencing of Table 1A before rounding indicated 
the presence of errors of about 2 in the tenth place—after larger errors had 
been removed. Accordingly, the error in Table 1A is about .52 in the eighth 
place. The error in Table 1B presents more of a problem. The larger interval 
made an appraisal by differencing considerably mo:e difficult; and, indeed, 
this difficulty accounts for our decision to present Table 1B to seven decimal 
places—instead of eight, as in Table 1A. Nevertheless, we think that errors 
larger than 2 in the ninth place would have been detectable, and that the table 
is therefore accurate to about .52 in the seventh place. 

Linear interpolation in Table 1A will produce better than five decimal accuracy 
throughout. In Table 1B, however, it will produce only three place accuracy, 
or a little better, at the lower end, improving to nearly five place accuracy 
at the upper end. To obtain full accuracy in either table will require the use of 
second differences, modified where necessary. Everett’s formula will be required 
from 7 = 1.4 to 7.0; elsewhere, Bessel’s formula will suffice. 


IIT. AN Example 


The distribution in Table 2 of hours to failure for 188 transmitting tubes, 
to which D. J. Davis (1952) attempted to fit an exponential curve, provides 
an instructive example of the use of our tables—illustrating the arithmetic 
procedures required as well as the awkwardness incident to grouped data. 
Ordinarily one would have no cause to complain about the grouping in Table 2: 
the groups are small, of equal length, and closed at both ends—at least on the 
arithmetic scale. On the geometric scale, however, the intervals are not equal, 
and the lowest one is wide open—extending all the way from 1.60954 = In 5 
to negative infinity. As the data in Table 2 stand, the geometric mean is little 
better than a guess. 

Although computing the geometric mean of a large number of ungrouped 
obsérvations is laborious, if performed manually, it is entirely feasible with 
modern computers; grouped observations, moreover, offer no facility for the 
geometric mean similar to the Charlier check for moments. When grouping 
appears unavoidable, we recommend: first, that the open-end interval be elimi- 
nated by recording the exact value of the smallest observation; second, that 
the class interval at the lower end be subdivided until no more than 1 percent 
of the population appears in the first cell. 

The lower part of Table 2 sketches the calculation of the arithmetic mean, 
the geometric mean, and 


y=nA — InG = 0.41276. 


Interpolation in Table 2 gives a value of yf = 0.55879 corresponding to y = 
0.41276; then 


é = 0.55879/0.41276 = 1.3538 





J. A. GREENWOOD AND D. DURAND 
TABLE 2 


Failure Data for 188 Transmitter Tubes 


Hours to Observed Hours to Observed 
Failure In X Freq. Failure X In X Freq. 


bq 


0-5 .91629 15 75-80 77.5 4.35028 


2.01490 18 80-85 82.5 4.41280 
2.52573 24 85-90 87.5 4.47164 
2.86220 30 90-95 92.5 4.52721 
3.11352 22 95-100 97.5 4.57985 
3.31419 15 100-105 102.5 4.62986 
3.48124 10 105-110 107.5 4.67749 
3.62434 110-115 112.5 4.72295 
3.74950 115-120 117.5 4.76644 
3.86073 120-125 122.5 4.80811 
3.96081 125-130 127.5 4.84812 
4.05178 130-135 132.5 4.88658 
4.13517 135-140 137.5 4.92362 
4.21213 140-145 142.5 4.95934 

: 4.28359 145-150 147.5 4.99383 


WMHOMNAEREWWHNR HE 
Anaananaaaananagnga&g;ca 
KOCOCOCONNKE RP WHEN NEE 


Zot 188 A = 32.20745 
> fin X = 575.174 In A = 3.47220 


sia 6055 | InG = 3.05944 
y 41276 


Source: Davis (1952) presented the distribution of failure times X in tabular form with a 
wide class interval and in graphical form with the above class intervals. The above 
table was derived from the histogram. 


and 


G@ = 32.20745/1.3537 = 23.79. 


We defer discussion of the standard error of 6 and 4 and of the question 
whether Davis’s exponential fit is acceptable to Section VII. 


IV. Tue ERLANG DISTRIBUTION 


The Erlang distribution, named after A. K. Erlang (cf. Brockmeyer, Halstrgm, 
and Jensen, 1948, pp. 25, 175) is a special form of the Type III with m = 0 
and integral values of p; it often appears in the literature of queueing theory. 
The process outlined above will provide maximum likelihood estimates for the 
Erlang distribution. If the calculated value of # is fairly close to an integer, 
one need simply accept this integer as the estimate of p and then substitute 
it in (3a)—with m = O—in order to obtain 4. If, however, the calculated value 
of # falls about halfway between two integers, one must substitute each integer 
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in (3a) and thence determine from (2) which pair of 6 and 4 gives the greater 
likelihood. 


V. THe GENERAL TyPkE III 


Although the procedure for estimating a and p for the gamma distribution 
can be extended to the estimation of a, p, and m for the general Type III, we 
cannot recommend it. Even for distributions where p is known to be greater 
than 1, we have devised no process that both exploits Table 1 and leads auto- 
matically to m. These likelihood estimates when found, moreover, are not 
sufficient—indeed no sufficient estimators exist—and maximum likelihood is 
therefore less attractive for three parameters than for two. If efficient estimators 
are satisfactory, one can well compute preliminary estimates by the method 
of moments and then improve these by using an algorism originally due to 
Fisher (1925), and extended to several parameters by Neyman (1949). 


VI. Tur Typr V DistrisuTion 


Since the Type V distribution does not command much interest at present, 
it suffices here to indicate how its parameters can be estimated. In analogy 
to (1), we write the Type V density 


ap (4g) eo [=47]. 


and find the log likelihood 
nPinA—ninI(P) —(P +1) } In(@a — M) — A > (e& — MD" 
and the estimating equations 


nP=A > (2 — MM)" (5a) 


-n< int) +nind = Dine - M) (5b) 


(P+1) OH (¢-M' =A DX (e@- MM” (5e) 


The application of Table 1 to the solution of (5a) and (5b) with M known 
should be obvious. 

The 3-parameter Type V may be treated analogously to the 3-parameter 
Type III; we again recommend the improvement algorism. In contrast to Type 
III, no difficulties arise when P < 1. 


VII. SranpArp Errors; Tests Or HyPoTHESEs 


Since 4 and # are sufficient estimators, their sampling distribution is in principle 
known (cf. Kendall, 1946, Sect. 17.31); but the marginal distribution of f and, 
a fortiori, the joint distribution of # and 4 are not suitable for computation. 
The second moments of the asymptotic joint normal distribution of d — a 
and # — p are obtained fairly easily by inverting a Hessian matrix. The process 
is implicit in Fisher (1922, p. 334-336) and is explicitly and conveniently given 
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by Kendall (1946, Sect. 17.46). From the matrix of the second differentials of 
the log likelihood 


Pp 1 
os | = a 
os # in Tp) 
dp* 
we find the Hessian determinant 


A= -4[1- Sinn | = — Fn 
a ? dp p a” 


and the variances and covariances 
n Var (6) = —1/n’ (6a) 
n Cov (4, p) = a/pn’ (6b) 
nm Var (4) = a'(1/p — 1/p'7’). (6c) 


Tables of n Var (@/a) and n Var () as functions of p appear in Masuyama and 
Kuroiwa (1951), for p as indicated in Section II. 

In large-sample work—indeed, formulas (6) are appropriate only in large- 
sample work—it is legitimate to replace a, p, and »’ by 4, #, and y’(); we may 
then evaluate the right-hand side of (6) by replacing 1/n’ by the first-difference 
approximation to 1/y’/—namely 


1/y’ = y' [yer — Yoro)/(yr — Yo) — A), 


where yo < y < y: , and y; — Yo = 0.01 (Table 1A) or 0.20 (Table 1B). This 
approximation is least accurate for 1.4 < y < 1.6; but even here it yields four 
significant figures in —1/y’. 

Equations (6) with a, p, and 7 not estimated from the data provide the equip- 
ment necessary for large-sample tests of relevant hypotheses. Since Davis repre- 
sented the data in Table 2 by an exponential, we may well wish to test the 
hypothesis p = 1 against the alternative p > 1. This test is a simple and in- 
structive special case of the problem of testing an Erlang distribution against 
the general gamma distribution. The test is performed by setting n = 188 
and p = 1 in (6a), yielding: 


Var (6 | p = 1) = 0.008282 
o(é | p = 1) = 0.0910 
p — p= 8.89c. 


Prima facie, we can reject the hypothesis p = 1 at the 0.001 level. 

Applying this test to the Davis example seems like stretching two points. 
First, as already mentioned, the error in the geometric mean due to grouping 
renders both f and its standard error dubious. Second, and much more basic, 
the figures in Table 2 are subject to a serious reporting bias, discussed by Davis, 
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that would tend to make them appear non-exponential even if they were essen- 
tially exponential. Accordingly, we must hedge the above significance level with 
suspicion. 

VIII. Potynomi1at AnD RatIonaL APPROXIMATIONS 


The approximations given in this section save space at the expense of com- 
putational convenience, and will probably be preferred by a statistician who 
is incorporating estimation of these parameters into a program to be stored in 
a computer. 


Maximum 
Error 


6 = y"'(0.5000876 + 0.1648852y — 0.0544274y’); 0 < y < 0.5772. 0.0088% 


, _ 8.898919 + 9.059950y + 0.9775373y? | 
6 = “9(17.79728 + 11.968477 y + y) 203772 Sy S170. —0.00547% 


Var 6 = p°(1.99629 — 1.21163y + 0.77255y*); 0 < y < 0.5772. 0.24% | 
Var 6 = £°(1.00513 + 0.89538 + 0.355584"); 0 < @ <1. 0.51% 
Cov (6, 4) = pa(1.99629 — 1.21163y + 0.77255y’);0 < y < 0.5572. 0.24% 


Cov (4, 4) = p4(1.00513 + 0.895386 — 0.35558f");0 < 6 <1. 0.51% 
Var @ = 4°(1.99780 + 0.73462y + 0.39306y*); 0 < y < 0.5772. 0.11% 
Var 4 = 4'(6"* + 1.00513 + 0.895386 — 0.355587);0<~<1. 0.22% 


These approximations were computed by the methods exploited by Hastings 
(1955). The error in # is less than o;/10 in samples of size less than 2.10°. 


IX, Accuracy REQUIRED 


One may ask when, if ever, the full accuracy of Table 1 is required for esti- 
mating some set of the parameters m, a, and p. Neither the quality of the data 
that one ordinarily encounters in practice nor the sampling errors that one 
might expect to incur in samples of reasonable size would seem to justify such 
accuracy. We find from (6a) that a sample of 10,000, which is seldom met in 
practice, yields a coefficient of variation in f# of 0.015 or more. Accordingly, 
calculation of p to three significant figures will ordinarily introduce a numerical 
error less than one tenth the sampling error, and this strikes us as adequate 
for confidence limits, decision functions, or tests of significance. 

One might, however, require considerably more accuracy in a problem calling 
for a complete job of curve fitting—including the estimation of theoretical 
frequencies in class intervals and the execution of some sort of test for goodness 
of fit. Owing to the substantial amount of calculation required in a curve-fitting 
problem, small errors introduced into the initial phases may grow into large 
errors toward the end. Thus, if one wishes to calculate x’ to four significant 
figures—and four figures are to be found in the common published tables of 
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percentage points—one cannot expect to get by with only four significant 
figures in the basic calculations of , 4, and #. Table 1 will provide just about 
the amount of numerical accuracy required to exploit the full accuracy of 
Karl Pearson’s table of the Incomplete Gamma Function, and it covers roughly 
the same range of p. 


X. Concitusion AND WARNING 


When m is determined from a source other than the data, the maximum 
likelihood estimators 4 and 6 are sufficient and have, therefore, a strong claim 
to being optimal. Moreover, with the aid of Table 1 or the approximation of 
Section VIII, the labor of calculating these estimates is not excessive—though 
certainly more than for moment estimates, since the geometric mean requires 
more labor than the standard deviation. 

For the general Type III there are no sufficient estimators. The labor of 
calculating efficient estimates, with the aid of the improvement algorism, is 
clearly greater than the labor of calculating moment estimates. It is not certain 
when the additional labor is worth the trouble. 

Finally, a word of caution is necessary. The methods described here yield 
maximum likelihood estimates only when the data constitute a sample drawn 
from a bona fide gamma population. When the data represent in fact a popula- 
tion o unknown character suspected of resembling the gamma approximately— 
and we fear this is the realistic interpretation—our estimators clearly are not 
maximum likelihood estimators: A and G estimate the population arithmetic 
and geometric means, but it is not clear that the calculated quantities 4 and / 
estimate any particular parameters. This does not mean that our methods are 


worthless for problems requiring graduation rather than estimation: quite the 
contrary, when data are to be graduated by a gamma distribution curve, our 
methods will, in general, provide a better fit by the Chi-square criterion than 
the method of moments—though not quite so good a fit as minimum Chi-square. 
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Experimental Designs to Adjust for Time Trends’ 
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Kingsport, Tennessee 


The problem of time trends is frequently present in chemical processes. Experi- 
mental designs which can be used in these situations have been developed by several 
workers. The designs of D. R. Cox are restricted to a few levels and can be used with 
both qualitative or quantitative variables. The randomly oriented designs of G. E. P. 
Box are useful when all of the variables are quantitative. Both of these types of 
designs presuppose the fitting of a first-order model to the data. Combinations 
of these designs can be made which allow the study of both qualitative and quanti- 
tative variables in the same experiment. Designs which permit the fitting of 
a polynomial response are available for one variable problem. Problems involving 
two variables can be handled by the proper combination of randomly oriented first- 
order designs. Satisfactory designs for situations involving more than two independent 
variables and a polynomial response of second order or higher are not presently avail- 


able. Further work in this area, will be of interest to experimenters in the physical 
sciences. 


INTRODUCTION 


Traditionally, the elimination of time trends has been accomplished by the 
use of blocking—Latin squares and randomized blocks—or by the use of co- 
variance analysis in which the covariate is a time trend. These techniques 
work satisfactorily in many cases. In the chemical industry, however, there are 
certain situations in which these devices are inadequate to provide all of the 
desired information. 

A situation which occurs frequently involves the study of heterogeneous 
catalytic reactions which exhibit a change in the activity of the catalyst with 
time. The chemist understandably would prefer to study the specific chemical 
reaction under conditions of constant activity of the catalyst. Thus, he is faced 
with the choice of replacing or regenerating the catalyst. Both of these pro- 
cedures are costly in time and money. In this situation, the randomly oriented 
designs of Box [1], which allow adjustment for time trends, are useful in pre- 
liminary investigations. 

In the application of heterogeneous catalytic reactions on a commercial 
scale, the operating conditions are frequently changed to maintain maximum 


1 Presented at the Conference on Statistics in Chemistry and Chemical Engineering, 
Gordon Research Conferences, New Hampton School, New Hampton, New Hampshire, 
August 24, 1959. 


2 Member of Northeast Tennessee Section of the A.S.Q.C. and the A.S.A. 
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productivity of the process. In this situation, it would be helpful to know if 
the effects of the controlled variables changed with time (time-variable inter- 
action). This information would be useful in setting the operating conditions 
and would be even more useful with the advent of computer control of chemical 
plants. ; 

In the discussion which follows, experimental designs proposed by G.E.P. Box 
[1], [2] and D. R. Cox [3], [4], [5], which allow adjustment for time trends will 
be presented. In addition, it will be shown that certain of these designs can be 
combined to permit the study of both qualitative and quantitative factors in 
the same experiment. For all the experimental designs described, the time 
intervals are equally spaced. 

It is hoped that this presentation will stimulate further work on other possible 
approaches to this interesting problem. 


Desiens or D. R. Cox 


Experimental designs which allow the study of the effects of controlled variables 
in the presence of time trends have been developed by Cox [3], [4], [5]. In these 
designs, the number of levels is very restricted (2, 3, or 4), but the variables 
may be either qualitative or quantitative. The trends usually are restricted 
to representation by a quadratic or cubic equation. 

As an example of a simple design, the plan for two treatments 7, and T, in 


TaBLzE 1 
Two treatments, four replicates 
(¢{=2,n =4,N = 8, p = 2) 
-5 -3 1 1 3 5 
7 1 -—3 —5 —5 -—3 1 
T: T; T: T: T T: T: 


&’ , &’ = values of levels of orthogonal polynomials 


four replicates is given in Table 1 where the notation (t, n, N, p) in this and 
subsequent tables gives 


t = treatments 
n = replications 
N = runs 
p = order of polynomial trend. 


In this particular design, the effect of the treatments is orthogonal to a quadratic 
time trend. Another similar design involving two treatments in six replicates 
is given in Table 2. These two designs have been designated by Cox as non- 
symmetrical. In a nonsymmetrical design, the arrangement of 7, and 7, when 
moving in the negative direction from zero along # (the orthogonal polynomial 
for the linear time effect) will be the reverse of that found when moving in the 
positive direction. In a symmetrical design, the arrangement will be the same 
in both directions. 
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TABLE 2 
Two treatments, six replicates 
(¢ = 2,n = 6, N = 12, p = 2) 


—7 —5 —3 —1 1 3 5 7 o. 
1 -17 -29 -35 -35 -29 -17 1 25 55 
T; T: T: T: Ti T: Ss. Fe. Te Ee 


Symmetrical designs which permit the estimation of treatment effects 
orthogonal to a cubic trend are available. For two treatments the smallest of 
these designs contains eight replicates. The positive half of this design is shown 
in Table 3. (The negative half of the design will be a mirror image of the half 


TABLE 3 
Two treatments, eight replicates 
(¢ = 2,n = 8, N = 16, p = 8) 


7 

-9 
—301 
qT; 


ll 

9 
—143 
qT; 


13 
21 
91 
TM 


15 
35 
455 
Ts 


shown.) In Table 4 are shown four designs of two treatments in eight replicates 
each. These designs are orthogonal to a quadratic time trend. It is interesting 
to note that these four designs considered together give a 2* factorial design 
where now each column vector gives either the two levels or indicates the 


TABLE 4 
Two treatments, eight replicates 
(¢ = 2,n = 8,N = 16, p = 2) 


&,’ 


35 
21 

9 
—1 
-9 
—15 
—19 
—21 
—21 
—19 
—15 
—9 
—l 
9 

21 
35 





70 HUBERT M. HILL 


presence or absence of four variables. The interactions of these four variables 
are mutually orthogonal, but orthogonal to only the linear component of the 
time trend. 

Cox also gives other designs which allowed the comparison of more treatments. 
These designs, however, have the treatments orthogonal to the odd-order 
polynomials only, while attempting to minimize the correlation between treat- 
ment effects and the quadratic component of the trend. Examples of several of 
the better designs are given in Tables 5 and 6. These are all symmetrical designs. 


TABLE 5 
Three treatments, six replicates 
(t = 3,n = 6,N = 18, p = 8) 


5 7 9 
—31 —22 —10 
—35 —42 —42 

T; Ts T2 
Ts T: T2 


TABLE 6 
Four treatments, four replicates 
({=4,n = 4,N = 16, p = 8) 


1 3 5 7 9 11 
—21 —19 —15 -9 -1 9 
—63 —179 —265 —301 —267 —143 

qT; T2 Ts Ts T. T; 


Analysis of these designs is straightforward in cases in which the treatment 
effects are orthogonal to all of the desired polynomials. In cases in which they 
are orthogonal to only the odd-order polynomials, it is necessary to solve the 
least squares equations to separate the treatment effects from the quadratic 
component of the trend. 


RANDOMLY ORIENTED DesieGns oF G. E. P. Box 


Experimental designs which permit the determination of the linear effects 
of quantitative factors independent of time trends of any order have been 
developed by Box [1]. In addition, it is also possible to re-estimate, up to a 
pre-determined order, the nature of the time trend. These designs are based on 
the mathematical model, given in Equation (1), which allows high order terms 
in the time trend, but only first-order terms in independent variables. 


(1) = Bo + Bik + Bath + <*> + Babh + ots + ante + +++ + ayty 


where ¢’ represents an orthogonal polynomial of order designated by the 
subscript. 


and x, represents the independent variables, i = 1, 2, -+- , p. 
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This model assumes that second-order terms and interactions among the p 
controlled variables are negligible and that the effects of these p variables do 
not vary with time (no interaction). 

The variation in the observations observed in an experimental study which 
involves a time trend, will be the result of the time trend, the controlled variables 
under study and errors assumed to be random. The variation due to the trend 
can be removed by fitting a polynomial to the data and the remaining variation 
studied for the effects due to the controlled variables. The designs of Box ac- 
complish this by choosing the levels of the variables x, , x2 , 73 , --* , %» So that 
the vectors in the design matrix are orthogonal to the vectors used to estimate 
the time trend. This is accomplished by choosing the design vectors so that 
they are linear combinations of orthogonal polynomials of higher order than 
those used for estimating the time trend. Further, since randomization is required 
to ensure valid inferences from data, the design vectors x, , X, , *** , X, are 
randomly oriented in the space defined by the higher order orthogonal poly- 
nomials. (This is necessary since the individual experiments comprising the 
completed design cannot be run in a random time sequence, time being one 
of the controlled variables under study). Finally, the p vectors are made mutually 
orthogonal. 

Examples of the construction of these designs are given in references [1] and 
[6]. The rest of this section will be devoted to a discussion of the results of 
possible alternative choices of higher order polynomials and of the combination 
of blocks to form new experimental designs. 


TABLE 7 TABLE 8 TABLE 9 


Ze 


Time Order 2 Xi Ze x1 Xe 


0.194 —0.749 —0.233 0.494 
0.959 1.189 —0.298 —1.113 
— 1.667 —0.744 1.830 0.735 
—0.975 1.054 —0.612 —1.040 
0.426 0.644 —1.363 1.117 
1.063 —1.394 —0.530 1.284 
1.063 —1.394 0.530 —1.284 
0.426 0.644 1.363 —1.117 
—0.975 1.054 0.612 1.040 
—1.667 —0.744 —1.830 —0.735 
0.959 1.189 0.298 1.113 
0.194 —0.749 0.223 —0.494 
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The design matrix in Table 7 is taken from reference [6]. The vectors in this 
design are orthogonal to a cubic trend and are linear combinations of poly- 
nomials of order 4 through 11 for N = 14. This design matrix will yield in- 
dependent minimum variance estimates of the coefficients in a first-order model 
(Equation 1). 

It is not necessary that the elements of the design vector be a linear combi- 
nation of all orthogonal polynomials higher than the time trend. If a linear 
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combination of the even-order polynomials is constructed, the design vector 
will be symmetrical in the sense described by D. R. Cox. This is illustrated by 
the design matrix in Table 8 in which the elements of the vectors are linear 
combinations of the polynomials of order 4, 6 and 8 for N = 12. On the other 
hand, a linear combination of the odd-order polynomials will result in vectors 
which are symmetrical with opposite signs. This is represented by the matrix 
in Table 9 in which the elements of the vectors are linear combinations of the 
polynomials of order 5, 7, and 9 for N = 12. 

It is interesting to note that in the design of Table 9 the estimates of the 
main effects of x, and x, will not be biased by any second-order terms in 2, 
and 2, . This is evident from the (X’X) matrix (Table 10) calculated for the 
following complete second-order model. 


(2) y= bo + dik + Dobe + Dats + ait, + Oe, 
+ O41} + Ay2%3 H+ Qy22%o + Oy:01k: + GreTok, + € 


Tasz 10 


a aun 
12 


2.9 


In Table 10, and in subsequent tables, only the diagonal elements and elements 
above the diagonal, of the symmetric X’X matrix are listed. 
The exact nature of the biases is best illustrated by constructing the alias 
matrix. Writing Equation (1) in matrix notation gives 
n= X,6; 
and for Equation (2) 
oo! X68; + X26. 


where X, represents the matrix of independent variables in Eq. (1), and X, the 
matrix of the additional independent variables found in Eq. (2) and where 
6, and 6, are the corresponding vectors of coefficients. If Eq. (2) is required to 
adequately represent the response function but only the terms of Eq. (1) are 
used, then the least squares estimates B, of §, will be biased and, as shown in 
Box and Wilson [7] 


(3) B, — 6, + [X/X,]"’X/X.6. or B,—>6, + AG 





EXPERIMENTAL DESIGNS TO ADJUST FOR TIME TRENDS 73 


where the arrow indicates that the quantity on the left is an unbiased estimate 
of the quantity on the right. The matrix A = [X/X,]~* X/X, is called the alias 
matrix. For example, using the results of Table 9 and 10, and assuming Eq. (1) 
is adequate when Eq. (2) is required, gives for the estimate of 6, 


(4) bs — Be + + [2.901 + 0.3022 — 1.7e12 + l.lan — 1.9e2] 
However, the full second order model in the 2’s is no longer fully independent 


of the time trend since the estimate of the second order polynomial constant 
in time is not orthogonal to the second order terms in z, and 2, . 


TaBLe 11 TABLE 12a TABLE 12b 


Block 2 Block 3 Block 4 
Time Order a(+)  2(+) a(—)  2(+) a(+)  2(—) 


0.223 —0.494 : —0.494 0.223 0.494 
1.113 : 1.113 0.298 —1.113 
—0.735 , —0.735 —1.830 0.735 
1.040 : 1.040 0.612 —1.040 
—1.117 F —1.117 1.363 1.117 
—1.284 ; —1.284 0.530 1.284 
1.284 , 1.284 —0.530 —1.284 
1.117 % 1.117 —1.363 —1.117 
—1.040 : —1.040 —0.612 1.040 
0.735 : 0.735 1.830 —0.735 
—1.113 : —1.113 —0.298 1.113 
0.494 ; 0.494 —0.223 —0.494 
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if a second block (Table 11) is run with opposite signs on 2, and z, , the 
interactions of linear time trend with the variables (z,t, and z2t,) are correlated 
with each other; but they are independent of the other terms in the model 
(Table 13). It should be recognized here that the time settings are in fact repli- 


TABLE 13 


bs a & au Pe) 
24 24 


5.8 0.6 
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cated, in each block, that is, if fresh catalyst were used at the beginning of the 
first block, fresh catalyst would again be used at the beginning of the second 


_ block. 


If two more blocks are constructed by changing the sign of only zx, or 2, 
(Table 12a and 12b), a design is constructed in which all main effects and inter- 
actions have zero covariances and only the quadratic terms are correlated (Table 
14). The precision of the estimates of the second-order effects will be related to 


TABLE 14 


bo b be bs a a ay 


48 48 
48 
11.6 


40.9 
34.2 
28.5 


the particular random orientation obtained in designing the experiment. Also 
it is probable that certain orientations will give lower correlation between 


the quadratic effects than others. The design in this example is obviously too 
large to be practical. Assuming a quadratic time trend and two variables, the 
smallest possible block would be six runs. Four blocks would then result in 24 
runs, which would not be an impossibly large experiment. 

In developing the four blocks in the previous design, the signs of the columns 
were changed in a factorial-type pattern (Table 15). This same device can be 


TABLE 15 


Configuration 


applied to larger experiments. With three factors, eight blocks would be required. 
With a minimum block size of six, 48 runs would be required. No reduction in 
the number of blocks is possible without confounding some of the desired effects. 
These designs are entirely too large above three variables, and some other 
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design should be devised to make it possible to use a quadratic model in the 
presence of a time trend. 


A ComMBINATION OF DeEsIGNns oF Box anp Cox 


In discussing the designs of Box and Cox, it has been shown that it is possible 
to construct designs which are symmetrical with the same signs or the opposite 
signs. It is evident that two vectors, one symmetrical with the same signs and 
the other symmetrical with opposite signs, must be orthogonal to each other. 
This can be illustrated by examination of tables of orthogonal polynomials. 
For instance in Table 16, which contains the orthogonal polynomials of order 


Tase 16 
&’ & 2" &3’ &1' 2! £2" fs" 
-3 +1 = ~3 cat" 
at -—1 +3 1 —3 


+1 —1 —3 —1 +3 
+3 +1 +1 +3 +1 


4 plus the é# and £££ vectors, the last two vectors have sums of zero. The 
orthogonality of ¢ with £ and éj is not due to the actual numbers present but 
to their arrangement. Therefore, changing the numbers in any of the first three 
columns will not affect the orthogonality of 1 and 3 with 2 as long as the required 


symmetry is maintained. Thus, we may combine the designs of Box and Cox 
to form new designs which will allow one to study both qualitative and quanti- 
tative factors in the presence of a time trend. As an example of this type of design, 
let us combine the design in Table 8 with the design from Table 2. This combi- 
nation yields the design in Table 17. 

The estimates of the effects of the factors will be orthogonal to a quadratic 
time trend (Table 18). Assuming a first-order model, the effects of x, and z, 


TaBLe 17 


Time Order 
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will be biased by second-order terms of z, and z, but not by any terms involving 
Z; . The effect of x, will be unbiased by second-order terms except for inter- 
actions of the linear time trend with z, and zx, . The interactions 2,2, and 22%; 
are uncorrelated with the quadratic trend but are correlated with the linear 
trend. Replication of this design by changing the signs of x, and 2, will give 
unbiased estimates of x, and xz, (Table 19). In addition, the interactions 2,275 
and 2,2; will be estimable independent of the linear or quadratic components 
of the time trend. 


Another design possibility is the reverse combination. An example with N = 16 
is given (Table 20). In this design, a linear model will yield biased estimates 


TABLE 20 
Time Order Ea 


—0.283 
—0.012 
1.100 
0.186 
—1.696 
—0.614 
1.493 
1.093 
—1.093 
—1.493 
0.614 
1.696 
—0.186 
—1.100 
0.012 
0.283 


SCONATh WN 


of the first-order effects. Replication of this design with the signs of xz; changed 
will eliminate these biases, and, in addition, will make z,x, and 2.23 estimable. 
These designs also can be used in situations in which no time trend is present. 


SEcoOND-OrDER Desiens IN ONE VARIABLE 


Another type of design, devised by Box and Hay [2], permits the estimation of a 
single factor response curve in the presence of a time trend. The model in this 
design is 

3. n = Bot Biti + Bote t--> $ Bits tard tad +--+ +a;d' 
where ¢ represents an orthogonal polynomial of order (7) and d represents the 
level of the factor under study. 

The problem in building a design to fit this model lies in picking the levels 
of the experiment in such a way that the a’s are always uncorrelated with the 


6’s. This requires that > td‘ = 0. The orthogonal polynomials of lower order 
are used to represent the time trend, and the levels of the factor (d) will be 
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a linear combination of higher order polynomials. Box and Hay have given an 
example of the construction and use of these designs [2]. 

In the situation in which the vectors are linear combinations of more than two 
higher order orthogonal polynomials, the equations become quite difficult to 
solve; but if a linear combination of only two higher order polynomials is used, 
they are readily solvable. Table 21 lists two solutions which are orthogonal to a 


TABLE 21 TABLE 22 TABLE 23 


3 


Zt 


Time Order d d’ d d’ 


—0.38 —0.95 
1.39 1.77 
—1.39 0.31 
—0.53 —0.87 
1.07 —0.98 
0.76 —0.39 
—0.76 0.39 
—1.07 0.98 
0.53 0.87 
1.39 —0.31 
—1.39 —1.77 
0.38 0.95 
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cubic trend. They were constructed as linear combinations of orthogonal poly- 
nomials of order 4 and 6 for N = 12. These vectors are symmetrical with the 
same signs. Table 22 lists two vectors which are linear combinations of poly- 
nomials of order 5 and 7 for N = 12. 

The quadratic component of the response is calculated as an orthogonal 
polynomial using the formula, 


3 
a, = @ — (24) a, -1. 


N 
In the design which is symmetrical with opposite signs, this reduces to 
dg =d; —1. 


It is interesting to note that in the designs in Table 21 the model can be ex- 
panded to include the terms ¢,d, . This interaction is orthogonal to t, , &% , di 
and d,. 

As originally used by Box and Hay, these designs were intended for the 
comparison of two response curves in the presence of a time trend. This was 
accomplished by running a block of two within each time period. The two 
variables were quantitative (dosage of a standard vs. an experimental chemical). 
This same device might also be used where one quantitative variable is being 
studied in the presence of a qualitative variable at two levels. 

This type of situation could also be studied by combining a symmetrical 
design with one of Cox’s designs as in Table 23. In this design, the effects of 
2, , X; , and 2, are free of biases (Table 24). The interaction 2,2, will be biased 
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TABLE 24 


bo b be a a au* 


12 
12 


ae 


* Calculated from d, = d,* — ()-d*/12)d; — 1. 


by the linear time trend. This bias may be removed by replicating the design 
with the signs of x, changed. When the design involving two blocks is used, the 
interaction t,7, and t,x, also can be estimated (Table 25). Designs of this type 
can also be constructed with quantitative variables having a design vector 


TABLE 25 


a 


- * Calculated from d,; = di? — (>.d3/12)d, — 1. 


which is symmetrical with opposite signs and a qualitative variable having a 
completely symmetrical configuration. However, in the smallest of such designs, 
N = 16. Replication of these designs with the opposite signs on 7, also allows 
the estimation of all the main effects and interactions. 


EXAMPLE 


This experiment was a study of a gas-phase reaction which was conducted 
over a solid catalyst. The reaction proceeded at elevated temperatures and was 
run continuously. The catalyst was activated during the first few hours of 
operation and the experimental study was started after the catalyst had leveled 
out. The “time” periods were measured in milliliters of reactant fed to the 
reactor. This system was used because the experimenters considered the decline 
in activity of the catalyst to be a function of the amount of material fed rather 
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than the time of operation. A uniform cut of material was allowed as a leveling 
out period between each set of reaction conditions. The product stream was 
analyzed by gas chromatography. Three responses were recorded: conversion, 
yield, and catalyst productivity. Only the results of the conversions will be 
shown. 

The design used was a randomly oriented design with four variables and 14 
sets of observations (see Table 26). This was constructed orthogonal to a fifth- 
order trend with the design vectors as linear combinations of polynomials of 
order 6-13 for N = 14. 

The design matrix with the resultant conversions is shown in the top half of 
Table 26. A second block with reverse signs was also run, and this is shown in 


TABLE 26 


Block 1 
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1 Missing value calculated by minimizing the residual sum of squares. 
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the bottom half of Table 26. The variables were coded so that the two-block 

experiment covered the range of conditions of interest to the experimenter. 

This caused the units recorded in Table 28 to appear to be rather unusual. 
The results of an analysis of the data are recorded in Table 27. The trend 


TABLE 27 
Analysis of Variance 


Regression Degrees of 
Coefficient Freedom Sum of Squares Mean Square 


—0.096 
1.710 
—1.216 
2.736 
% —0.358 
Residual 
Total 


— 
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could be easily represented by a linear equation. There appeared to be a difference 
between the two blocks, probably due to a difference in initial activity of the 
catalyst attained in the reactiv: ion step. Variable x; showed the largest effect 
followed by z, and x, . Variable z, also showed a small but significant effect. 
Accordingly, a path of steepest ascent was calculated (Table 28). Runs were 
made at points 29 and 30. The conversions obtained were 24.4 and 26.2%, 


TABLE 28 
Path of Steepest Ascent 
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respectively. These compare favorably with the calculated values of 24.8 and 
28.6%. The results obtained from the other responses mitigated against a further 
move along this path, so this portion of the problem was dropped and further 
work was concentrated on other aspects of the problem. 

I am indebted to Dr. E. L. McDaniel and Dr. H. 8. Young for the experimental 
results shown. 
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Tests for the Validity of the Assumption that the 
Underlying Distribution of Life is Exponential.’ 


Part I 


BENJAMIN EpsTEIN 
Wayne State University and Stanford University 


It is frequently useful to test, on the basis of life test data, whether or not one 
is justified in assuming that the underlying distribution of life is exponential. This 
paper, which appears in two parts, describes a number of graphical and analytical 
procedures for testing this assumption. Part I of the paper contains descriptions of 
the mathematical and graphical procedures. Part IT contains several worked examples. 


1. A GRAPHICAL PROCEDURE 


Let us first note a graphical procedure, which is particularly useful if one 
has a large amount of data. Suppose that the underlying distribution is really 
exponential, then the cumulative distribution function (c.d.f.) F(é) is given by 


(1) F(t) = 0, t<0 
=1-—-e, #t>0. 
It is further assumed that @ > 0. From (1), it follows that 


(2) y = log (=) = t/6. 


That is, if we plot y = log [1/(1 — F(é))] against é, we get a straight line with 
slope 1/@. This suggests the following graphical procedure for testing departures 
from exponentiality: Suppose we place n items on life test and let t; < 4 < --- 
< 4, be the times when the items fail. In plotting the graph, the value of F(é;) 
associated with ¢,; will be 


. 


a 

This is done because E[F(t;)] = i/(n + 1). If the exponential assumption holds, 
then the plotted points can be fitted well by a straight line passing through the 
origin. If experimentation is discontinued at time ¢, (r < n), when the rth failure 
occurs, then we expect a good linear fit up to ¢, . 


1 Research sponsored in part by the Office of Naval Research. 
83 
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Remark 1: It may happen that the underlying density function of life is the 
two-parameter exponential p.d.f. 


® Kts0, y= te", t>A>d0 


= 0, elsewhere. 
The associated cumulative distribution function is 
(5) F(t; @, A) = 0, t<A 
=1—e“%4" t*>A>0. 
From (5), it follows that 


(6) log (Sires! = (t — A)/6. 


Hence, if the underlying distribution is two-parameter exponential, we will 
still get straight lines if we plot log [1/(1 — F(¢,))] against ¢; . However, instead 
of passing through the origin, the straight line cuts the ¢-axis at the point t = A. 


Remark 2: Suppose that one starts the life test at time ¢ = 0 with N, items. 
Let N(t) be the number of items surviving after a time ¢ has elapsed. Then, if 
the exponential assumption holds, N(é) = Noe~‘”’ and log N(#) will plot linearly 
against ¢. This is an equivalent and alternative graphical procedure to the 
one given above. 


2. THE x° TEST FOR GOODNEss OF Fir 


If one has an extensive body of failure data, then one can estimate the param- 
eter 6, divide the time axis into a number of intervals, use the estimate of @ to 
find the expected number of failures in each interval, and then use a x’ goodness 
of fit test with (k — 1) degrees of freedom, where k is the number of intervals 
(or classes) into which the time axis has been divided. More precisely, let 6 be 
the best estimate of @ based on n failures and let the time axis be divided into 
k intervals by the (& — 1) times 4, < & < +--+ < &-1 . Furthermore, let 0; be 
the observed number of failures in the zth interval and let e; be the expected 
mumber of failures in the ith interval, where e; = np; and where 


ta ‘+ 
@ p= [ger ae, 
0 O@ 


” 1 -2/6 ° 
p= | 7 dz, for #=1,2,---,k-—1 
bina 


Pr — 1 o-s/ dz. 


In terms of the above definitions, x’ is defined in the usual way as 





TESTS FOR VALIDITY. | 


_ . (0; -— e)". 

(8) x = a . 
The statistic on the right-hand side of (8) is distributed as x’ (k — 1) (chi-square 
with (&K — 1) degrees of freedom) if the sample size (number of failures) is large. 

It is well known, however, that the chi-square goodness of fit test has several 
drawbacks. Among them are its large sample character and dependence upon 
the choice of the number and position of the intervals into which the time axis 
is divided. For a discussion of this point, see Gumbel [20] and [21, p. 28]. Ac- 
cordingly, it seems worthwhile to give statistical procedures which are com- 
paratively free of this kind of difficulty, and this is done in the sequel. 


Remark: For an excellent survey of the chi-square goodness of fit test, one 
should read Cochran [12]. 


3. A CRITERION BASED ON THE CONDITIONAL DISTRIBUTION 
or Torat LIvEs 


This test makes use of basic properties of Poisson processes (for a discussion 
of Poisson processes, see Feller [18]). It can be shown (see Appendix 1) that if 
one observes a Poisson process for a fixed length of time 7 and if r events occur 
in [0, T] at times 7, < 72 < «++ < 1, < T, then these times (after being sub- 
jected to a random permutation) can be considered as r independent observations 
on a random variable uniformly distributed over (0, 7). For r even moderately | 
large, >-3., 7; is approximately normally distributed with mean r7’/2 and 
variance r7”/12. This can be used to test whether or not the data are drawn 
from a Poisson process. 

One can also show (see Appendix 1) that if one observes a Poisson process 
until exactly r events occur (where r is a preassigned integer) and if the events 
occur at 7, < t2 < -** <17,, then the (r — 1) random variables 7, , 7. , --- , 
T,-, can be considered, when unordered, as (r — 1) independent observations 
on a random variable which is uniformly distributed over (0, 7,). For r moderately 
large, >-1z! 7; is approximately normally distributed with mean (r — 1)r,/2 
and variance (r — 1)r?/12. As before, this fact can be used to test whether or 
not the underlying distribution is Poisson. 

In the context of life testing, it is clear that the comments just made for 
Poisson processes go over directly in case items on test are replaced immediately 
by new items. This gives rise, of course, to a Poisson process with constant rate 
and we are b. - to the situation just treated. If it happens that we are dealing 
with a non-replacement situation (where failed items are not replaced), then 
all we need to do is to use total lives T(r;), rather than the failure times, 7; . If 
T(r;) is the total life observed in getting the ith failure, then T(r:) < T(r.) < 
+++ < T(r,). In the case where the life test starts with n items and is terminated 
at a preassigned total life 7*, the number of failures observed will be a random 
variable r with 


T(r,) =n ,T(m) = 1 +@—1)n,-::, 
Tir) =ntrates+tr1tM@—rti)r,. 
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As before one can show that the total lives T(7,), T(r2), «++ , T(7-) can be 
considered as being drawn from a density function which is uniform over [0, 7*]. 
If the life test ends as soon as the first r failures occur, then the (r — 1) random 
variables 7'(7;), T(r2), --* , T'(7--1) can be considered as being drawn from a 
density function which is uniform over [0, 7'(7,)]. 

The fact that the conditional distribution of total lives is uniform over suitable 
intervals makes it quite evident that one has a good tool for detecting whether 
the failure rate is indeed constant (as it must be in the pure exponential case). 
Thus, for example, the contamination of a purely exponential distribution by 
early failures would manifest itself in a pronounced tendency to get too many 
failures clustering together in the early part of the time or total life axis, thus 
violating uniformity. If the underlying distribution is really described by a two- 
parameter exponential with A > 0, then we should expect to get too few failures 
occurring in the early part of the uniform distribution, thus violating uniformity. 
If the failure rate changes, for example increases with time, then this should 
result in a tendency for failures to cluster together as time goes on, again violating 
uniformity. For a case in point, suppose that half-way through the life tes+ 
the failure rate increases by a factor of two; then one can expect to find roughly 
twice as many points in the second half of the interval as in the first and this 
would violate uniformity. 

If the amount of failure data observed is quite small, then we can expect to 
detect only fairly large changes from exponentiality by the preceding tests. If 
we have a substantial amount of data, one can use a x’ test to detect whether 
the conditional distribution of times to failure or total lives (whichever is 
appropriate) deviates excessively from being uuiform. Some illustrations of 
how this is done could be found in the examples given in Part IT. 

We now give some tests which make essential use of the fact that total lives 
between successive failures are-mutually independent and drawn from a common 
exponential distribution. 


4. A Test ror ABNORMALLY EARLY FAILURES 


It may happen that the underlying failure distribution is exponential. How- 
ever, the first one or two failures may for one reason or another occur abnormally 
soon. How do we test for this kind of deviation from the exponential? 

Suppose that 7, < 72 < +++ < 7, are the first r failure times and that we wish 
to test whether 7, is abnormally small. To carry out the test we recall [15], [17] 
that if all the 7; are drawn from a common exponential, then 7'(7,), the total 
life in [0, 7,], and 7'(7, — 7;), the total life in [r, , r,], are distributed independently 
of each other. 27(7,)/0 is distributed as x*(2) and 2T(r, — 7;)/@ is distributed 
as x°(2r — 2). Therefore, (r — 1)T(1:)/T (7, — 1:1) is distributed as F(2, 2r — 2). 
Clearly we may reasonably assert that 7, is abnormally small if the ratio is too 


small. More precisely, if a is the significance level of our test, we will say that 
7, is abnormally small if 


(9 T(r) < Fee 
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where F* is a lower a@ point of the F(2, 2r — 2) distribution, i.e., 


(10) Pr(F(2, 2r — 2) < FE) =a. 


In a similar way we can obtain a criterion for detecting whether both 7, and r, 
are abnormally small. As before, we know that if all the 7; are drawn from a 
common exponential distribution, then T(7.), the total life in [0, 7.], and 
T(r, — 72), the total life in [7, , r,], are distributed independently ofeach other. 
Furthermore, 27'(7,)/@ is distributed as x’(4) and 27(r, — 12)/@ is distributed 
as x’*(2r — 4). Therefore, the ratio (r — 2)T(r2)/2T(r, — 72) is distributed as 
F(4, 2r — 4). If the ratio is too small, we may assert that both 7, and 7, are 
abnormally small. More precisely, if a is the significance level of our test, we 
will say that both 7, and 7, are abnormally small if 


(11) T(r) < Pet PE. =) 


where F** is a lower @ point of the F(4, 2r — 4) distribution, i.e., 


(12) Pr(F(4, 2r — 4) < F&*) =a. 


5. A Test ror AN ABNORMALLY LONG First FAILurE 


If it happens that the underlying distribution is really two-parameter ex- 
ponential (see the density (4)) with minimum life A > 0, then we get a deviation 
from the exponential which can be detected by virtue of the fact that the time 
when the first failure occurs is abnormally long. A reasonable test for this kind 
of departure from the one-parameter exponential is to work with the ratio 
(r — 1)T(r:)/T(r, — 11). This statistic, which compares the total life (suitably 
weighted) in (0, 7,] with the total life in [7, , 7,], is distributed as F(2, 2r — 2) 
if A = 0. If we wish to test the hypothesis that A = 0 against the alternative 
A > 0, then we wiil reject A = 0 at the significance level a, if we get the abnorm- 
ally long first failure time 7, or if 


(13) T(r) > os He 


where G* is an upper a point of the F(2, 2r — 2) distribution, i.e., 


(14) Pr[F(2, 2r — 2) > Gt] = a. 


6. A Test FoR WHETHER THE MEAN LIFE (or RatTE oF FAILURE) 
IN THE First HAtrF or A Lire Test Dirrers SIGNIFICANTLY 
FROM THE MEAN LiFe (or RaTE OF FAILURE) IN THE 
Seconp Harr or A Lire TEst 


This is a rough and ready kind of procedure which can be used to detect 
gross changes over time. Comparing the first half of a life test with the second 
half of a life test is intended only to indicate the kind of procedure that one 
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may use. The actual choice of which two periods to compare is often dictated 
by a@ priori considerations. 

Suppose that 7, < tz. < +++ < 72, are the first 2r failure times. Let us compare 
what happens in [0, 7,] with what happens in [r, , 72]. If the underlying distri- 
bution remains exponential with a fixed @ throughout, then 7(7,), the total 
life in [0, 7,], and T'(r2, — 7,), the total life in [r, , r2,], are distributed independent- 
ly of each other and 27(r,)/@ and 27(r2, — 1,)/@ are each distributed as x°(2r). 
Therefore, 7'(7r,)/T(t2e — 7,) is distributed as F(2r, 2r) and we can use the F 
table to give us tests of significance. Whether these tests should be one-or two- 
sided depends on what kind of change we are trying to test for. 

More generally if we have the first s failure times r, < 7, < +++ <1, , we 
may wish to compare the period [0, 7,] with [7, , r,] where k < s. In this case 
T(7,), the total life in (0, 7,) and T(r, — 7,), the total life in [7, , 7,], are dis- 
tributed independently of each other. If the observations are drawn from a 
common exponential distribution with fixed 0, then 27(r,)/@ and 27(r, — 7,)/@ 
are distributed as  x°(2k) and x’(2s — 2k), respectively, and therefore 
(s — k)T(r,)/kT(r, — 7) is distributed as F(2k, 2s — 2k). The F table can 
then be used to give tests of significance. Results in [16] are relevant to the 
above considerations. 


7. A Test ror WHETHER OR Not THE MEAN LIFE FLUCTUATES 
DURING THE Lire TEST 


In procedure (6) we tested whether or not the mean life (or failure rate) in 
the first half of a life test differs from that in the second half. A natural extension 
is to consider _ one has kr ordered failure times 71; < t12 < *-** < Tir < Ta < 

+ < top S++ S te < +++ < MH, arranged in k groups each of size r. Assuming 
that the mean life 6 is sida within each group, the object is to test whether 
or not @ fluctuates from group to group. 

It is easy to show that one can apply Bartlett’s test for homogeneity of variance 

[6], [7]. The connection comes from the fact that T(7,,), T(ta — Tir), - 
T'(tne — Ta- 1)e)y the total lives in (0, Tir]; [n1, ’ Tor], en Tee [ru- -1)r» Teel; respectively, 
are mutually independent and 2[7'(7;, — t¢;-1)-)]/@ is for each j = 1, 2, ,k 
distributed as x’(2r). Furthermore, T(m,) = T(tir) + T (tor — as + +. + 
T (Tie — Te-1)r) is such that 27'(r,,)/@ is distributed as x’(2rk). It is then very easy 
to verify that 


2ri{ log Tee) T(r) 7 [log T(71,) + log T( 72, — tir) Fees + log T(t — Ta-1,)] 
k al 1 


+ “ae 


(16) * an mem 


‘€ 
The hypothesis that @ is the same over the entire group is rejected if x’(k — 1) 
is too large. If x2(k — 1) is the upper a point of x’(k — 1), ie., Pr(x'(k — 1) > 
x2(k — 1)) = a, then if the left-hand side of (16) exceeds x 2k — 1), we will 
aa the hypothesis that @ is the same throughout, at the diapdticies level a. 
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8. A Specrat Cask or ProcepureE (7) WHEN r = 1 


Suppose that in the preceding test described by (16) we let r = 1, and let 
T11 = 71,721 = 72, °** , aNd tT, = 7, thus getting 7, < 72 < ++ < %. T(r), 
T(t. — 71), «++ , T(t, — Te-1) are mutually independent and 27(r; — 1;-:)/0 
is for each j = 1, 2, --- , & distributed as x°(2). Further, T(r.) = T(r) + 
T(t. — 11) + +++ + T(t, — 7-1) is such that 27(r,)/60 is distributed as x*(2k). 
If we make r = 1 in (16) we get 


2I{ log 2) me i [log T(71) + log T(t2 — 71) +--+ + log T(n — nd} 

14 e+) k -t 1 

(17) ~ - — 1). 
Remark 1: The numerator of (17) can be written as 


“a » log {T(r; — r;-1)/T(t; — 15-1}, 


T(t; — 4-1) = 2 T(t; — 1;-1)/k = T(11)/k. 


T(r; — 75-1), j = 1, 2, +++ , k are total lives between failures and T(r; — 7;-_,) 
is the arithmetic mean of these values. This test can be viewed as checking 
whether or not @ is constant from one observation to the next. It is interesting 
to note [27] that the likelihood ratio test obtained when one tests the hypothesis 
H, , that the distribution of total lives between successive failures is described 
by the p.d.f. f(x; 6) = [1/6] e-*”*, against the class of alternatives H, , that the 
p.d.f. is given by 


f(x; 0, B) "nae * eae 


(note that if 8 = 1, f(z; 0, 1) = f(x; 6) = [1/6] e~*”’), leads to (17) as the test. 
Again one is led to the rejection of the hypothesis that the underlying distri- 
bution is exponential with constant @ throughout time if the left-hand side of 
(17) becomes large. The test based on (17) can be viewed as an alternative 
to test 2, which is based on the chi-square test of goodness of fit, or to test 3, 
‘which is based on the uniformity of the conditional distribution of ordered 
total lives. In this connection, we would like to return for a moment to test 3. 
It is well known that if an observation u is taken at random on a variable dis- 
tributed uniformly in the interval [0, 1], then — 2 log u is distributed as x*(2). 
From this it follows immediately that 


~2 ¥ leg [Bed | 


t=1 
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is distributed as x’(2r) where 7* is a preassigned total life and the number of 
failures, r, is a random variable. In a similar way it can be seen that if one termi- 
nates life testing after a fixed number, r, of failures, then 


Si x ~ Re 2 


is distributed as x*(2r — 2). These tests are similar in nature to Bartlett’s test. 


Remark 2: The reader may wonder whether he should use formula 17 or the 
procedure just described to test whether the life test data deviates from the 
exponential distribution. Actually it depends on which class of alternatives 
one is testing against. Thus formula (17) (which is based on total lives between 
failures) appears te be best when testing whether the distribution of total lives 
between successive failures is exponential with constant @ against the alternative 
that the p.d.f. is Type III (see Remark 1) or is a member of the Weibull family 
i.e., of the form 


fla; 0, y) = (ya Oe", =2>0,0>0,7>0. 


Note that if y = 1, f(x} 0, 1) = f(z, 6) = [1/6]e*”". 

If, however, we view the sequence of failures and hence associated total lives 
as being generated by a stochastic process, then the assumption that total lives 
between successive failures are drawn independently and at random from an 
exponential p.d.f. with constant mean life @ is equivalent to assuming that the 
observations are being generated by a Poisson process with constant rate \ = 1/8. 
Bartholomew has shown in a recent paper [4] that if one tests the hypothesis 
that a stochastic process is Poisson with constant rate \ against the alternative 
that it is time dependent Poisson with rate \(é) = A(At)* (note that when a = 0, 
A(t) = for all ¢), then one is led to procedures inyolving actual times of events 
and not times between events. In the context of life testing we are led to tests 
of the form — 2 >>1z? log [T(r,)/T(r,)], i.e., tests involving total lives and not 
total lives between failures. 


9. A Test BASED ON THE Maximum F DistrisuTIon 


A quick test for homogeneity is based on some results of Hartley [22]. Consider, 
e.g., the situation in test 7 where we have kr time ordered failures which we 
arrange in k groups each of size r. Then a test for homogeneity is obtained by 
computing 


_. max T(11,) T (Te, => iage:.©* * a i = T h~1)*) : 
(18) U min {T(71,), T (tae — Tir) °° 9 T(tee — Tee-12)} 


Clearly one should reject the hypothesis of homogeneity if U is too large. Hartley 
has given 5% points for U for various values of k and r and has made it plausible 
that the test is almost as efficient as the Bartlett homogeneity test. For the case 
where r = 1, we simply have k failure times 7, < 72 < +--+ < 7, and U becomes 


(19) — max aT T1 T(t. — Ti)>5 T Ter ~~ Tr- 
~ min {T(71), T(12 — os oo+ SU(m — T-1)} 
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This is a quick alternative procedure to that given in test 8. Exact 5% and 1% 
points have been worked out for this case by Hartley [22]. 


10. Tests FoR ABNORMALLY LONG PERIODS IN WHICH 
THERE Is NO FAILURE 


Test 5 gives a reasonable procedure for detecting whether or not it takes too 
long for the first failure to occur. Now we mention a test for whether any of 
the intervals between successive failures is too long. In connection with some 
work on tests of significance in harmonic analysis, R. A. Fisher [19] has obtained 
the distribution of a very useful statistic. Stating his result in the terminology 
of life testing, let us define the statistic Z as 


(20) Z=a max {T(7), T(t. ne pe See T(t. ae Tn-1) } 
T (12) 


where 7, < Tz < +++ < 7, are the first n ordered failure times. It has been shown 
by Fisher that 


r 
= 


(21) Pr(Z > Z.) = X ("ya — kZ,)"", 


where r is the largest integer less than 1/Z, . The hypothesis of homogeneity 
is violated if Z becomes too large. Fisher gives values of Z, such that 
Pr(Z > Zo) < a for a = .05, .01 for values of n up to 50. 

(21) is the distribution of 


n 
max {x.}/ >> Vis 
l<isn i=1 


where 2, , 2, *** , Z, are m independent observations made on a random variable 
X having p.d.f. [1/é]e~*”’, x > 0, 8 > 0. In the life testing problem we identify 
: 2 with 27 (r,)/0, where T(7,) is the total life in the interval [0, 7,]; x. with 
27 (r2 — 71)/0, where T(r, — 7;) is the total life in [7, , 72]; x; with 27 (7; — 7-1) /0 
and x, with 27(r, — 7,-:)/6. Thus 


n 


max {x;}/ >. X; 


1sisn t=1 
becomes (20). 


Another test for abnormally long intervals between failures, which includes 
Fisher’s test as a special case, follows from a paper by Cochran [11] dealing with 
the largest of a set of estimated variances. A small table of 5% significance 
levels in Cochran’s paper was extended in a paper by Eisenhart and Solomon 
[14] which appeared as Chapter 15 of Techniques of Statistical Analysis. In 
addition, a table of 1% significance levels is given in [14]. 

Cochran’s test deals with the following: We are sampling from an underlying 
normal distribution with variance o*. We are given k independent variance 
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estimates of o”, 8}, --- , 82, each based on » degrees of freedom. Cochran defines 
the statistic g; by the equation 


k 
(22) 9 = 8/208, 4=1,2,-+-,k. 


In [11] he gives the distribution of max; {g;}. 
Since »s?/o* is distributed as x*(v), Cochran’s result includes, as a special 
case, the distribution of 


k 
max {2,}/ >) 2; 
1sisk #=1 


where 2; , 22, °** , % are k independent observations made on a random variable 
distributed as x’(2r). In the life testing problem we identify x, with 27(r,)/@, 
where 7'(7,,) is the total life in the interval [0, 7,,]; x. with 27 (72, — 11-)/0, 
where T (Ter =— Tir) is the total life in [rT ’ Tarl; zx; with QT (ris -_ T-1)r)/8, 
and XT with 2T (rir are T¢e-1)r)/0. Thus 


k 
max {z,}/ >) 2; 
1sisk t=1 
becomes max; {g;}, where g, is given by 
(23) gs = Trip — te-ry)/T(t), t= 1,2, +++, ke 


In using the tables in [14] we will act as if we have k independent estimates of 
variance, each of which is based on » = 2r degrees of freedom. The test based 
on (20) is a special case of (23) with k = n and » = 2. 


11. A GrapHicaL ProcepurE BASED ON THE 
KOLMOGOROV-SMIRNOV TEST 


Consider a life test which is terminated after a preassigned total life 7*. The 
number of observed failures will be a random variable r and let these failures 
occur at times 7, < 72 < «++ < 7, with associated total lives T(7:) < T(r.) < 
+++ < T(z,). We have seen previously that 7'(7,), T(72), «+> , T'(r-) are (when 
subjected to a random permutation) uniformly distributed over (0, 7%), if 
the underlying distribution is exponential with constant 0. Following Barnard 
[1] and Cox [13], a simple graphical way of checking this hypothesis is to plot 
a random walk diagram, in which the proportion of failures (relative to r) that 
occur at or before 7'(r) is plotted against 7'(r). If the hypothesis being tested 
is true, then the random walk diagram should fluctuate around the straight 
line joining the points (0, 0) and (7*, 1). Using the results of Kolmogorov [24] 
and Smirnov [31], one can base a test of significance on the size of r? D, , where 
D, = maximum difference (measured vertically) between the plot of the random 
walk and the straight line. 


Remark: Analytically the situation is as follows. We can define the sample 
c.d.f. F,[T(r)] as 


(24) rF,[T(r)] = number of failures occurring at or before total life T(r). 
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The theoretical c.d.f. is 


(25) F[T(r)] = T(7)/T*, for 0 < T(r) < T*. 
Thus if D, is the Kolmogorov-Smirnov statistic, then 
(26) rD,=r max | F,(T(r)) — FiT(>)]| 


0<T(r)<T* 
rT'(r) 
T+ | 


If r is large (> 80), the hypothesis that the underlying distribution is exponential 
with constant 6 will be rejected on the .10 level if the maximum deviation as 
measured by r! D, exceeds 1.22. The .05 and .01 rejection values of r* D, are 
1.36 and 1.63, respectively. For small values of r one should use the tables 
computed by Z. W. Birnbaum [10] and Massey [26]. 

In the preceding, we considered life test data which arise from life tests which 
are terminated after a preassigned total life 7*. Similarly, suppose that we carry 
on a life test until (n + 1) failures occur at times 7, < r2 < +++ < Tas, . Then 
we know that 7(r,), T'(r2), «++ , 7'(7,) are uniformly distributed over (0, T'(7,+1)), 
if the underlying distribution is exponential with constant @. Hence it is obvious 
that a graphical procedure analogous to the one just given involves replacing 
r by n and T* by T(r,,1). The test of significance can be based on n' D, , where 
D, = maximum difference (measured vertically) between the plot of the random 
walk and the straight line joining the points (0,0) and (T(r,.:), 1). 


= max _ | (number of failures up to T(r)) — 
0<T(1r)<T* ‘ 


Remark: Analytically one defines F,[T(r)] as nF,[T(r)] = number of failures 
occurring at or before total life T(r). The theoretical c.d.f. is 


F[T(r)] = T(t)/T(ta+1.), for 0 < T(r) < T(ta41). 
Thus if D, is the Kolmogorov-Smirnov statistic, then 
nD,=n max  |F,[T(r)] — F(T(7)) | 


O<T(r)<T (ta 41) 


- mex (number of failures up to T(7)) — as) 


O<T(r) <T(ra41) T(Ta+1) 


In the foregoing we have discussed two-sided tests. There are frequently 
situations where one should use a one-sided Kolmogorov-Smirnov type test. For 
a discussion of such tests and approximate tables, see Birnbaum and Tingey 
[9]. Numerical examples illustrating both two-sided and one-sided procedures 
will be given at the end of the article. 


12. A Test BASED ON CONDITIONAL RATE oF FAILURE 


A characteristic property of the exponential distribution is that the con- 
ditional probability of failing in the interval (¢, ¢ + A ¢) given that it has survived 
up to time ¢ is independent of ¢. This is immediately verified in the case where 
f(t; 0) = [1/o]e*”, t > 0,0 > Oand Pr(T < t) = F(t; 0) = 1 — e”. For this 
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case, the conditional probability of failing in (¢, ¢ + A #) given survival in (0, t) 
is given by 


f(t) At _ 
ik 1— FO ~ 


The practical significance of this remark is that if we start with a large number 
N, of items on life test, divide the axis into intervals (0, 7’), (7, 27), (27, 37’), -+- 
(where T' is suitably chosen) and if m, , m2 , m3 , -** are the numbers of items 
failing in these intervals, then 


Se cehac ai ciepeaneetilia 5s”. <iaagenietnehnabaiaisan ee ctnaliinitains 
N’N-—-1,’N—-1—7n,’ 7N—™ — NM — °°? — Mp1 


should fluctuate within reasonable limits about a constant value, namely the 
failure rate. 

In case the underlying distribution is not exponential, f(t)/[1 — F(é] will 
be some function g(t) and solving the equation 


(28) f()/[1 — FO] = 9d) 
and recalling that f(i) = dF(é)/dt, one finds that 
(29) Fi) =1-—e°™ 
where G(t) = fo g(r)dr. 


Remark: If G(i) = t/@ (i.e., g(t) = 1/6), one gets the one-parameter exponential 
case. If G(t) = 0,¢ < A and Gif) = (t — A)/0,t > A (ie, gD) = 0,1 < A, 
g(t) = 1/0,4 > A), one gets the two-parameter exponential. In case it happens 
that G(é) is not a linear function of ¢, then the sample conditional probability 
distribution will not fluctuate about a constant, but will instead approximate 
the function g(t). Thus, for example, if G(#) = at", the sample conditional 
probability distribution should approximate a curve of the form Ct*~’ (or a 
straight line if one plots 


MH ee I a Se ee 
{log , log 3" — log 7—™* —, \ 


against 
{ LT jog oe. log oe sh 
5 ? og 2 ? og 2 ’ 


Some GENERAL COMMENTS 


A dozen procedures have been given for testing for departures from an ex- 
ponential distribution of life. Some of these procedures are essentially equivalent 
whereas some procedures are superior to others in certain circumstances. In 
what follows we discuss various aspects of these tests which may be helpful 
to the reader. 
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Tests 1 and 2 have already been discussed adequately in the main body of 
the text. An alternative to test number 1 is test number 12. Substantial sample 
sizes (of at least 50) are needed before we can expect to detect substantial 
departures from an exponential distribution of life, by means of tests 1, 2, and 
12. Tests 3 and 11 are somewhat more sensitive than 1, 2, and 12 and will be 
more effective in detecting changes from exponentiality. Also we might expect 
to use these tests even for sample sizes as low as 10, although it should be recog- 
nized that for small sample sizes the power of the test will be low. 

Tests 3 and 11 both stem from the same fundamental properties of Poisson 
processes. It is difficult to compare the power of test number 3 with that of 
test number 11, since the power will depend on the set of alternatives against 
which one is testing the null hypothesis (that the underlying distribution is 
exponential with constant mean life). Test number 11 does have the advantage 
that it can be carried out graphically. There is one situation where test number 
3 can be shown to be “best.”’ This happens when one is testing the hypothesis 
that successive failures are being generated by a Poisson process with constant 
rate \ against the alternative that it is time-dependent Poisson with rate \(é) = 
ae’*. Cox [13] has shown that under these circumstances test procedure 3 is 
best. In [3] Bartholomew compares the process of 3 and 11 when the alternative 
is of the kind just described. In view of [13], procedure 3 is, of course, more 
powerful then procedure 11, but the difference in power is small. 

Test procedures 7 and 8 are direct consequences of Bartlett’s test for the 
homogeneity of variances. It has been shown by Moran [27] that the test pro- 
cedure 8 (formula 17) is an asymptotically most powerful test of the hypothesis 
H, , that total lives between successive failures are distributed according to the 
p.d.f. f(x; 0) = [1/é]e-*”", > 0, 6 > O, against the alternative that they are 
distributed according to a Type III distribution of the form 


e as Se —1 —2/8 
f(z; 6, 6) me” é . 


Bartholomew [5] has carried out calculations which compare the power of test 
procedure 8 (Bartlett’s test) with that of two other test procedures against 
various kinds of alternatives. He finds that test procedure 8 is the most powerful 
of the three procedures studied against alternatives which are either Type 
III or Weibull. 

Tests 7 and 8 make essential use of total lives between events, whereas many 
of the other tests involve the total lives themselves. That tests based on total 
lives should play an important role is intuitively clear because of the fact that 
the conditional distribution of total lives observed up to a preassigned total 
life 7* or up to a preassigned number of failures is uniform. The precise form 
taken by the tests depends on the alternatives against which one is testing. 
The interested reader should see papers by Bartholomew [2], [3], [4], [5], and 
Cox [13] on this point. 

In [2] Maguire, Pearson, and Wynn present data on time intervals between 
mine accidents. In [2] and subsequent commentaries on it [1], [4], the data are 
subjected to a variety of tests for departures from an exponential distribution 
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with constant mean. These papers give additional insight into the relative power 
of various procedures. 


SUMMARY 


In this paper we give a variety of procedures for testing, on the basis of life 
test data, whether there are significant departures from an exponential distri- 
bution of life. The particular procedures that one should adopt depend on the 
class of alternatives one is testing against. A number of the tests are based in 
an essential way on fundamental properties of Poisson processes. Questions 
involving choice of tests are considered, and a number of examples are worked 
out. 


APPENDIX 1 


(a) We wish to prove the assertion that if one observes a Poisson process 
for a fixed length of time T and if r events occur in (0, T) at times 7, < 7. < --- 
<7, , then these times can be considered, when unordered, as 7 independent 
observations on a random variable uniformly distributed in (0, 7). 


Proof: Let t; = time when the 7th event occurs; then 


Pr(r, < 4 < 1m t+ An ,72 <b < tm + An,°°:, 


t <t, <7, +Ar,|r events 


occur in (0, 7')) 


r r —-\T . st 
= ne? TI Ar, QMe _ ot! PAs, (4 Sane Se. 


i=1 r! oe T" i=1 

Hence the times {7;} considered as unordered are uniformly distributed in 
(0, 7). An equivalent statement is that the set {7,/T},7 = 1,2, --- , ris uniform 
over (0, 1). 

(b) We wish to prove the assertion that if one observes a Poisson process 
until exactly r events occur (where r is a preassigned integer) and if the events 
occur at 7; < tT. < +++ <7,, then the (r — 1) random variables 7, , 72, +++ , T=1 
can be considered, when unordered, as (r — 1) independent observations on a 
random variable which is uniformly distributed over (0, 7,). 


Proof: Let t; = time of occurrence of the 7th event, then 
Pr(r < t <1+ An, <& <2 +An,°°: ? 


Tr-1 = tp-1 < Tr-1 + At,-1 | Tr < l, mn Tr + Ar,) 


= Xe" Il Ar; Qe) 6-arry Ar 
— + oe . 


er ‘253 
ais a i ne 4 2 4 Ss 
r i=1 


T —_ 
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Hence the times 7; , 72 , --* , Tr-1 considered as unordered are uniformly 
distributed over (0, 7,). An equivalent statement is that the set {r;/7,}, 7 = 
1, 2, --» ,r — 1, is uniformly distributed over (0, 1). 

In life testing one can immediately apply the preceding provided that one 
replaces the times 7, by total lives T'(r,). 

The results that we have proved are used in an essential way in test procedure 
3 (a criterion based on the conditional distribution of total lives) and in test 
procedure 11 (a graphical procedure based on the Kolmogorov-Smirnov test). 


APPENDIX 2 


We should like to make more specific the connections between formulae 
(16) and (17) (tests 7 and 8) and work by Neyman and Pearson [29], Nayer 
[28], and Bartlett [6], [7] on the homogeneity of k variances. 

The L, test devised by Neyman and Pearson [29] is designed to test the 
following: We are given a sample of n; observations from each of k normal 
populations. We should like to test the hypothesis H, that all k populations 
have the same variance, but may have different unspecified means py; , 7 = 1, 2, 
--» , k. If the ith sample mean is given by 


a4 ni 
Oi De 2i,5/N 
7-3 


and the 7th sample variance by 


b= ¥ (es — 20m; 


i=1 


if, further, we define s2 by the formula 


Ah 
s = > n3i/N, 


t=1 


where N = )-*_, n, , then it can be shown that the appropriate test to use is 


‘ based on 
ne 


i=1 \Sy 


Rejection of the hypothesis occurs when Ay, becomes sufficiently small. 

The distribution of Ay, is difficult to find, but it is possible for large n; to 
work with (Az,)””. For details on the form of the limiting distribution one 
should refer to M. G. Kendall [23, pp. 298-299]. 

For the special case where all the n, are equal to a common n, (Az,)””” reduces 
to the L, statistic of Neyman and Pearson. In this case (Az,)””” becomes 


k ¢ 1/k 
t= 11(%)", 
t=1 ‘a 


since N = kn. Nayer has given 1% and 5% values of L, for n = 2(1)5(5)20(10)50 
and k = 2, 3, 4, 5, 10, 20, 25, 50 in [28]. 
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What is the connection between the L, test and tests 7 and 8? In test 7 one 
has kr time ordered failures 1: < t12 < *** S tie S Ta S + S Tar Se 
< Th < +++ < Te arranged in k groups each of size r. If 6 is constant throughout, 
then T (tis), T (Tor = Tir), Be ery T (Tre = Tn-1)); the total lives in 
{0, rie], [71+ » Tar], -** » [Ta-1ye » Ter], are mutually independent and 2[7(7;, — 
T-1)r)]/@ is for each j = 1, 2, --- , k distributed as x*(2r). But this means 
that our problem is equivalent to a situation where samples of size n = 2r + 1 
are drawn from each of k normal populations, and we wish to test for equality 
of variances (but where the populations may have different unspecified means 
Ms). The LT, test is thus obtained by replacing 8 by T (ti: — T4-1)e) and replacing 
< by T(r, ees T¢s-1)9) where 


a k 
T(r — T-1)r) = 2 T (tir ae: T 4-1) r)/k — T(1:)/k. 


Thus L, becomes 


T (Te) /k 


Computationally it is easier to take logarithms and so to find L, we would 
first compute 


i. es I [2s iar cucu) | 


ee T(t, 
log L, = 7 2» log T( tir — Teu-1).) — log (inw)). 


Remark: In entering the L, table one should take the sample size as n = 2r + 1 
and the number of samples as k (the number of groups). 

While the L, test, as we have defined it, is suitable for most practical purposes, 
a somewhat better test, in the sense that it is unbiased, is to work with the 
analogue of Bartlett’s test. Bartlett [7] has shown that if we define y”” as 


k 2 vi/y 
vs; 
ge Te 
vm > v8 


i=1 


where »,; = n; — Landy = )-*_, », , then (— 2 log u)/c, where 


citer pth C— 3) 


is approximated well by the x7(k — 1) distribution. In the life testing context, 
n, =n = 2r+1,», = 2r for each i and » = )-*_, », = 2rk. Hence c becomes 


1 k 1 
cit aga a 


Further, »””” = yw/"* = L, . Hence (— 2 log »)/e becomes (— 2rk log L,]/c, 
which is approximated well by the x(k — 1) distribution. But this is precisely 
formula (16). 

In the case where r = 1, we are dealing with k failure times, 7, < 7. < -:: 
< 7, , and one is concerned with the k mutually independent random variables 
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T(r); T(t. - 71); Pe T(t, —_ Tr-1)5 where 2T (r, —_ T;~1)/0 is for each 7 = fe 
2, -++ , k, distributed as x*(2). L, becomes 


Il im 7..." sa 
sei ( T(1)/k 
where 


T(r.) = T(1:) os T(t _ 7) oh eee - T(t aie Te-1) and further =i ee Es 
6k 

is approximated well by the x*(k — 1) distribution. This is precisely formula 

(17). 


Remark: If the L, table is used, take the sample size n = 3 and k equal to the 
number of failures. Rejection of the hypothesis that @ is a constant throughout 
the life test occurs for large values of x*(k — 1). 

Hartley’s maximum F test [2] forms the basis of our test number 9. This 
test can be looked upon as a “quick and dirty” analogue of test 7 and 8. 

Suppose that one is given k normal populations, having unknown variances 
o, ,i = 1, 2, --- , k, and that we wish to test the hypothesis that o? = 03 

- = of. To do this we are given variance estimates s 7, each based on » degrees 
of freedom, for each oj . Then Hartley’s test involves computing 


Fax = Max 8;/min 8; 


and rejection of the hypothesis (that all of are equal to a common (but 
unknown) o”) occurs if Faz is too large. In the life testing context we have 


the following obvious analogue. Suppose we are given k total lives {7;}, 7 = 
1, 2, --+ , k, each of which is based on r failures and suppose that we wish to 
test whether @ is constant throughout. Then Hartley’s test becomes F... = 
max 7',/min 7; and one enters his tables with number of samples equal to k 
and with the degrees of freedom for each sample, v = 2r. 


APPENDIX 3 


If one observes a Poisson process with constant rate d for a fixed time 7 and 
if r events occur in (0, 7’) at times 7; < 72 < +--+ < 1, < T, then the random 
variables u,; = 7,/T, i = 1, 2, --- , r can be considered, when unordered, as r 
independent observations on a random variable uniformly distributed over 
(0, 1). In a similar way if one observes a Poisson process until exactly r events 
occur (where r is preassigned) and if the events occur at rt, < t2 < +--+ <17,, 
then the (r — 1) random variables »,; = 1;/r, ,4 = 1, 2, «++ ,r — 1 can be con- 
sidered as (r — 1) independent observations on a random variable which is 
uniformly distributed over (0, 1). Now the {u,;} or {v,;} can be considered as 
arising from random division of the unit interval. Hence it is clear that work 
by P. A. P. Moran [27] and by Barton and David [8] on problems connected 
with random division of the unit interval are relevant to our present consider- 
ations. As is usual, if we are in life testing situations we will work with total 
lives rather than failure times. 
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AppENDIX 4 


It is known that if t, , 4, --+-+ , ¢, are intervals between successive occurrences 
of a random process, then the statistic 


o= >> |t; —i|/Qni, where i= Dot; 
t=1 i=1 
can be used to test the hypothesis that the underlying process is Poisson with 
constant parameter \. This statistic has been studied by Sherman [30] and 
recently it has been applied by Bartholomew [2], [5] to test for departures from 
an exponential distribution. 

In the life testing situation mw can be computed as follows. Suppose that we 
observe a life test until » events occur. Let 7, < 72 < --- < 7, be the times 
when the n events occur and let the associated total lives be T'(71) < T(12) < 
+++ < T(r,). The total lives between failures are T(7,), T(72 — 71), T(73 — 72), 
+++, T(t, — T,a-1). Thus w becomes 

oe PG. = 4g oe 


i=1 n 


/2T(r,). 


It would appear from Bartholomew’s work that although a is not most powerful 
against certain specific alternatives, it is quite good against a broad class of 
alternatives. This is very important since some of the procedures we have 
discussed are “best” against certain alternatives, but quite poor against other 
classes of alternatives. Thus it would seem that if we know a priori that the 
departures from exponentiality are likely to be of a certain kind, then we may 
want to choose the test procedure which is best against the particular class of 


alternatives in question. If, as is often the case, such a priori information is 
lacking, then w is a good all purpose test. 
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Programming Fisher’s Exact Method of 


Comparing Two Percentages 


W. H. Rosertson 
Sandia Corporation 


Exact solutions to many statistical problems are now possible through the use of 
high speed computing equipment. This paper describes the application of a high 
speed computer for determining the exact probability statement associated with the 
problem of comparing two percentages. 


INTRODUCTION 


When two percentages are to be compared to see if a significant difference 
exists, there are several statistical methods for making the comparison. One 
method which is not used to its full advantage is Fisher’s Exact Method (1). 
Approximate methods are more frequently used because the computations re- 
quired by Fisher’s method are generally time-consuming, even when a desk 
calculator is available. 

With the advent of digital computers, Fisher’s method should be employed 
more often. In fact, now that small compuiers are becoming generally available, 
the determination of exact probability statements through the use of computers 
is a concept that the entire statistical world should be considering. 

In this short paper, a method of programming Fisher’s method is given. This 
program has been used on the Royal McBee LGP-30* computer and can be 
adapted easily for use on many other computers. 


FisHEer’s Exacr Meruop 


To illustrate the use of Fisher’s method, assume that the two percentages 
to be compared are in the fraction form p, = a/(a + 6) and p, = c/(c + d) 
where a = min (a, b, c, d). Here, a, b, ; c, d are the success-failure dichotomies 
representing two classifications which we wish to compare. This may be more 
easily illustrated in a tabular form: 


Outcome 
Classification Total 
1 a b a+b 
2 c d c+d 
Total a+e b+d a+b+c+d 


* The LGP-30 is a desk-size digital electronic computer having a 4096 word magnetic 
drum memory and utilizing a command structure with 16 commands. Except for subroutines, 
the subject program utilizes 188 storage locations. 
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The exact probability of observing the two fractions p, and p, when there is 
no class difference, i.e., the true dichotomy for each class is in the proportion 
(a + c): (b+ d), is given by 


p, = @tyierO'!at+orb+a! (1) 
. albic!d!(a+b+c+d)! 


To obtain the final probability to use in judging whether there is a significant 
difference, we must add to (1) the probabilities of more divergent fractions 
than those observed. The next more divergent situation (assuming p, < p, 
is obtained by decreasing a and d, and increasing b and c by unity. Thus, the 
two-way table for the second term of the series is: 


Outcome 
Classification Success Failure Total 


1 @m&®=a-—1 b&b =b+1 a+b 
2 G@ =c+1 d=d—1 c+d 


a+e b+d a+b+c+d 


Se (a+b!e+a!a+o!(b+_9)! (2) 
*~ @-DYIG+YV!E+V)V!d-VY!i@+b+e4+ oO! 
In general, 
P.= a@+bd'lc+AO!a+olb+d)! (3) 
*"@-t+)b+i-VYie+ti-VYid-i+)la+b+ce4+_d! 


fort = 1,2,---,a+1. 
Finally, the probability of observing two fractions as much or more divergent 


than p, and p, when there is no difference in the sampled populations is given 
by 


a+l 


P=) P,. (4) 


#=1 
EXAMPLE 


As an example of the application of Fisher’s Exact Method, consider a situa- 
tion where a silicon dip coating for vacuum tube protection is to be evaluated 
in terms of the difference between observed failure rates of the protected and 
the unprotected tubes. A total of 690 tubes of the same type and quality were 
subjected to the same environment. The following success-failure dichotomies 
were observed: 


Outcome 


Failures Non-Failures 


Protected 2 338 
Unprotected 7 343 


Total 9 681 
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The failure rate observed for the protected tubes (p, = 2/340 = .0059) was 
less than that observed for the unprotected tubes (p. = 7/350 = .0200). The 
experimenter wishes to determine if the difference between these failure rates 
is more than could be attributed to chance, thereby indicating that the pro- 
tected tubes are superior. The decision was made to judge significance at the 
10% probability level. 


By equation (1) the probability of obtaining the exact results observed when 
there is no difference is given by 


340! 350! 9! 681! 


P. = 513381 7! 343! 690! 


Using logarithms of factorials the computations for this result are as follows: 


log 340! 714.7076 
log 350! 740.0920 
log 9! 5.5598 
log 681! = 1635.4344 

3095.7938 
log 2! 0.3010 
log 338! 709.6460 
log 7! 3.7024 
log 343! 722.3097 


log 690! = 1660.9612 
3096.9203 


Log P, (iifference) = 8.8735-10 
Log P, = .0747 


Similarly, the probability of observing failure rates where the protected tubes 


are even more superior to the unprotected tubes, are given by equation (3). 
Thus, 


_ _ 340! 350! 9! 681! 
~ 11339! 8!342! 690! 


_ _ 340! 350! 9! 681! 
0! 340! 9! 341! 690! 


P, .0189 


P; .0021 


By equation (4) 
P = .0747 + .0189 + .0021 = .0957 


and it may be concluded that the failure rate for protected tubes is just barely 
significantly less than that for unprotected tubes. 
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THe ProGRAM 
Equation (1) may be reduced to a form which is more conducive to program- 
ming. 
P,=| a+b  _a+b-1 
* Latb+t+etdat+b+c+da-—1 
at+b+c+da-—2 b+c+d+1 
[oteote—lote-2 c+1) 
a a-—l a-—-2 1 
{ ct+td ct+d-1 c+td-2 d+i | (5) 
b+e+db+c+d-—-1lb+c+d-—2 b+d+1 


or 


_F (a+b) —2z Wa@to-y F_C€+9 -2z 

nies Hetitet+s ~2 I a-y Nete+ a -- © 
This lends itself to a three loop program, where the bracketed terms in (5) 
are a self modifying loop. The calculation of P, , --- , P, is a repetition of 
the program with the elements incremented as in (2). It may be verified that 
P.41 (a; = 0) is equivalent to the last bracketed term in (5). Thus, to compute 
P.,; , only the third loop of the program is required. Further, if the program 
does not provide for skipping to the third loop when a; = 0, the computer will 
attempt to divide by zero. A flow chart for the program is shown in Figure 1. 


CALL IN DATA; a, c, (a+b), (c+) 


COMPUTE AND STORE: b, d, o-y*@, otb-x = a+b, ate-y = ate, 
ct+d-z* ctd, b+ctd-z=btctd, atbtctd-z = a+bicté 


COMPUTE (a+b-2/(e+b+c+d-) 


[_ COMPUTE (a+e-9)/(0-y) | 
[ MT BY A a STORE 
[TEST eye | 
No [ves ] 
DECREMENT o+b-x BY | | DECREMENT otb-y BY | | 
DECREMENT otb+c+d-x BY | 


ADO Pj TO EP; 
DECREMENT o, TO a; —I 
INCREMENT ¢, TO ¢ +1 [ NO} TEST FOR NEGATIVE YES | PRINT =P, | | STOP } 
Figure 1—Program flow chart for comparing two percentages by Fisher’s exact method. 
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CoMMENTS 


In the program outlined, the accumulator may be required to operate with 
a large range of numbers, the largest of which is 


@+bd!iat+ol(b+c+ad! 
al blec!(a+b+c+d)! 


Since rounding error can be a consideration (due to the large number of repeated 
division-multiplication operations) scaling may be a problem if the program 
is written in fixed point. In floating point, rounding error should not affect the 
results since the final probability is seldom required beyond three decimal places. 

The computing time required for this program varies directly with the size 
of a and c but is independent of the size of b and d. That is, the speed of the 
program depends largely on the number of division-multiplication cycles in- 
volved, and the total number of these is (a + 1) (3a + 2c)/2. 

The program could easily be modified to stop when (4) reached a preset 
significance level. Thus, in most comparisons where the fractions are not signifi- 
cantly different, the computing time would be reduced. 
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Misclassified Data from a Binomial Population* 
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The method of maximum likelihood is employed to estimate the binomial pa- 
rameter p whenever an erroneous observation or report results in recording c defective 
samples when there are actually c + 1. The proportion of such erroneous observa- 
tions is also estimated, and asymptotic variances and covariances of the estimates 
are obtained. A numerical example is included. 


1. INTRODUCTION 


Inspection procedures cften involve the selection and examination of a single 
sample consisting of » individual units. If the number of defective units in a 
sample exceeds a specified acceptance number c, the production lot represented 
is unacceptable. If the number of defectives is c or less, the lot is acceptable. - 
Due to a natural reluctance on the part of some inspectors to reject a lot because 
of “only one too many defectives” in a sample, samples containing c + 1 defec- 
tives are sometimes reported as having only c defectives with the result that 
lots are accepted which would have been rejected had prescribed inspection 
procedures been rigidly followed. When this situation exists, any estimate of 
the overall average percentage (or fraction) of defective items in production 
presented for inspection which fails to properly take these misclassifications 
into account will be in error. 

In this paper, we employ the method of maximum likelihood to estimate the 
process or population fraction p when samples are subject to the type of mis- 
classification described above. Here, we are concerned with a binomial population, 
but the methods employed are quite similar to those used previously in studying 
related problems of misclassification in samples from Poisson populations [1], 
[2}. 

Suppose the number of defectives actually present in a sample of size n is a 
binomial random variable with parameter p. Let 6 designate the probability 
that a sample which actually contains c + 1 defectives is misclassified by report- 
ing it as containing only c defectives. In all other cases we assume correct 
observation and reporting of defectives. Let X designate the number of defec- 
tives reported in a sample of size n. The probability function of random variable 
X, f(x) = Pr(X = x), may be written as 


* Sponsored by the Office of Ordnance Research, U. S. Army. 
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(ra a = 0,1,---(—1),¢€+ 2),::: 


ta {yu w+ Ce2}. a=0 


a n e+1 a n-e-1 we 
(1 a, 4p a-pr, zac, 
where 0 < p < 1,0 < 0 < landec = 0,1,2,---,n—1. 


2. DERIVATIONS 


Consider N observations of X. Let r, be the frequency with which X = zx 
occurs. It follows that >-* r, = N, and the likelihood function of these N obser- 
vations may be written as 


P(t, , +++ tw 5D, 8) 


= 1+ QE) a -oe T1Y{Cra-ar}". @ 


Taking logarithms of (2), differentiating with respect to p and @ in turn, and 
equating to zero yields the estimating equations 


r.(n — c)@ 


eh 
+G=plet 0d — 9 Fe om 


aL r(n — c)p ( Tear ) i, (3) 


a6. ¢€+DA—-pD+m—om \i-@ 


where L is written for ln P. . 

The required estimates # and 6, when they exist, will be found by simul- 
taneously solving these two equations. We follow the customary notation in 
this paper of employing (*) to distinguish estimates from parameters. 

The pair of equations (3) is equivalent to 


Oresi/(l — 6) + > xr, = Nnp, (4) 
rn — c)/[(e + 11 — p) + (nm — Op] = r.4:/p(1 — 8). 


We solve the second of these equations for 6, substitute into the first, and 
simplify to obtain 


pnN|n — 2c — 1] + p\nN('+ 1) —-(n — 2c — 1) ; xr, —1,(n — c) 


— Terie + »} —(¢+ n> a, — ren} = 0, (5) 
6 = {r. — terile + IL — p)/(n — ep} /. + te41)- 
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The first of these equations is quadratic in p, and the required estimate # 
must then be a positive real root in the closed interval (0, 1). With # thus de- 
termined, 6 follows from the second equation of (5). 
When no misclassification has occurred, and when therefore @ = 0, The first 
equation of (3) yields the well-known result 


p= > ar,/N. 


In the event that (7) all sample observations are c’s, in which case r, = N’ 
Tex: = 0, and >—* zr, = Ne, or (ii) all sample observations are (¢ + 1)’s, in 
which case r, = 0, rou: = N, and >-* ar, = N(c + 1); then maximum likeli- 
hood estimators # and 6 will not exist, since in each of these cases a simultaneous 
solution fails to exist for the two equations of (3). In the event that r, = r.,, = 0, 
no estimate is possible for @, but from the first equation of (3), we would have 
f = > ar./nN. Although of theoretical interest, these exceptions seem 
unlikely to be important in practical applications. 

It is possible to construct other samples for which equations (5) fail to yield 
acceptable estimates of p and 6. However, for large N such samples will be 
very improbable and their occurrence in practical applications should be in- 
terpreted as a suggestion that perhaps probability function (1) is not applicable 
to the random variable actually observed. 


3. SAMPLING Errors oF EsTIMATES 


The asymptotic variance-covariance matrix of (, 6) is obtained by inverting 
the matrix consisting of elements that are negatives of expected values of the 
second partial derivatives of Z. Differentiating (3) and taking expected values, 
we obtain 


E(—#L/dap’)/N a eae (i ra Jae fe {1 


B(—0'L/ap 80)/N = on = — ED = on = E(-0'L/00 ap)/N, ©) 


se Ss 
E( °L/a0)/N = se" eA? G ak 6)?’ 


where A = [(c + 1) (1 — p)/(n — c)6 + pl, f. = JO), ferr = fle + 1), and 
915 = 94i(D, 8) with 7, j = 1, 2. 
The asymptotic variances and covariance follow as 
V@)~ = ¢22/N i922 — $12); 
V() ~ ¢u/NGu¢2 — Giz); (7) 
Cov (p, 6) ~ —¢1:2/N@u¢2 — Gia). 


When desired, the coefficient of correlation between estimates may be computed 
as 


13.8 ~ —P12/V Pirgaa- (8) 
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The variances and covariance given by (7) and the correlation coefficient given 
by (8) are applicable in all cases where maximum likelihood estimators # and 6 
exist. 


4. Spectat CAsE IN WuHicH c = 0 


In this case, estimating equations (5) become 
Nn(n — 1)p’ + {nw —(n—1) Dar, — nn — rhe _ {> ar, — nh =, (9) 
0 0 


6 = {ro — r,(1 — p)/np}/(ro + 1), 
and the ¢;; of (6) become 


| ~~» @ >i {i no MI 

me — p 1—p+npo \lp—1— p+npe 
_ - s  S 

Pi2 (1 —» + mpl ’ 


Goo = np(l — py | eB + 7 aaa | 


[1 — p + npey 


When 6 = 0, we obtain the well-known result for the variance of #: 
Vp) = pl — p)/mN. 


5. AN ILLUSTRATIVE EXAMPLE 


To illustrate the practical application of the foregoing results, we consider a 
hypothetical example consisting of the reported number of defectives in 1000 
samples of 40 items each. For purposes of this illustration we assume that 

= 2. Data for this example are tabulated below. 


Number Defectives Frequency 
z Ts 


296 
366 
273 
34 
25 
5 

1 


Summarizing, we have N = 1000, n = 40, c = 2, pe ar, = 1145, r, = 273, 
and 7.4, = 34. On substituting these values into equations (5) and solving 
we obtain 


f = 0.0299 and 6 =0.61. 


These estimates are to be compared with population values p = 0.03 and 
6 = 0.605 which were used in constructing this example. 
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In order to calculate estimate variances and covariance, we first employ 
(6) to calculate 


G11 = 13.39, Giz = —1.94, and Y22 = 0.59, 


where p and @ have been replaced by their estimates 0.0299 and 0.61 respectively. 
From (7) and (8) it follows that 


V(p) = 0.000 000 75, V(6) = 0.0017, 
Cov (f, 6) = 0.000 002 5, and 5 = 0.07. 
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Book Reviews 


Norsert L. Enricx, Quality Conirol—A Manual of Quality Control Procedure 
Based Upon Scientific Principles and Simplified for Practical Application in Various 
Types of Manufacturing Plants. 3rd ed. (New York: The Industrial Press, 1959), 
184 pp. Price $4.50. 


This is a revision of a book that first appeared in 1948 and was revised and enlarged in 
1954. The second edition added a part II containing additional materiel on control charts 
and “analysis of variance’ and the third edition made minor corrections and revisions and 
enlarged the chapter on “analysis of variance”. 

The book is a very simple account of the principles of statistical quality control intended 
primarily for the practical man. It presents a series of simple procedures that could easily be 
followed by the uninitiated. The discussion is largely in terms that could have immediate 
meaning to managers and engineers concerned with quality control and the procedures are 
justified largely by appeal to general knowledge and common sense rather than mathematical 
proof. 

The chapters run ‘Fundamentals of Inspection” (4 pp.), “Procedure in Installing Lot- 
by-lot Inspection” (18 pp.), “Sampling Continuous Products” (3 pp.), “Installing Process 
Inspection” (11 pp.), ‘Special Control Charts for Use when Equipment is Old and Worn” 
(5 pp.), etc. Although the book is intended for beginners, those already versed in statistical 
quality control may be interested in some of the special techniques used to simplify its presen- 
tation. Statisticians working in the industrial field will find the chapters on Tolerances and 
Allowances and Mass Production Gaging helpful in enlarging their understanding of the 
industrial phases of these subjects. 

Although the extreme simplicity of the book will be a help to those first seeking light on 
modern quality control, they should be warned that the subject is really much more sophisti- 
cated than presented. They should also be warned that the simplicity introduces inaccuracies 
and possibly even misunderstanding. It is hoped that readers will not be induced to install 
- a system of quality control without first consulting an expert. 

A feature of the book that the reviewer found particularly annoying was the use of the 
term ‘Allowable Percent Defective”. In the last chapter this is explicitly identified with the 
Acceptable Quality Level, but when first presented it appears to have the characteristics 
of the AOQL. On p. 13, for example, we read “... the Allowable Percent Defective is only 2. 
Thus, occasionally individual lots having between 2 and 3.3 percent defectives will be passed 
through (not considering those bad lots which slip by now and then due to the “luck of the 
draw” errors normally present in the sampling process) although in the long run the overall) 
proportion of defectives will be held at 2 percent.” (Italics the reviewer’s) If the Allowable Percent 
Defective is merely the AQL, this statement is not correct. 

The author also offers acceptance sampling plans for process control and control charts 
for acceptance sampling without clearly distinguishing the basic nature of the two tools. 
His acknowledgement that many of the techniques presented assume normality of output 
comes rather late in the book and his “analysis of Variance” is more “components of variance 
analysis” based on ranges than what is usually called analysis of variance. But these are 
thoughts that disturb the sophisticated reader. It is doubtful that they would bother the 
beginner. If the latter is led by the author’s simplified approach to dig further into the subject 
to enlarge his knowledge and correct his initial misunderstandings, the book will have fulfilled 
its purpose. A, J. Duncan 


115 





















TECHNOMETRICS 





Fesruary, 1960 


NOTICES 


SouTHERN REGIONAL GRADUATE SUMMER SESSION IN STATISTICS 
AT UNIVERSITY OF FLORIDA 


The 1960 Southern Regional Graduate Summer Session in Statistics will be 
held at the University of Florida at Gainesville from June 20 to July 29, 1960. 
The University of Florida, North Carolina State College, Virginia Polytechnic 
Institute and Oklahoma State University have agreed to operate a continuing 
program of graduate summer sessions in statisics to be held at each institution 
in rotation. The first such session was held at Virginia Polytechnic Institute 
in the summer of 1954. 

It is the purpose of this program to serve: (1) teachers of introductory statisti- 
cal courses and college teachers of mathematics who want formal training in 
modern statistics; (2) research and professional workers who want intensive 
instruction in basic statistical concepts and modern statistical methodology; 
(3) professional statisticians who wish to keep informed about advanced special- 
ized theory and methods; (4) prospective candidates for graduate degrees in 
statistics; and (5) graduate students in other fields who desire supporting work 
in statistics. 

The session will last six weeks and courses will carry three semester hours of 
credit. Not more than two courses may be taken for credit at any one session. 
The summer work in statistics may be applied as residence credit at any one 
of the cooperating institutions, as well as certain other universities, in partial 
fulfillment of the requirements for a graduate degree. The program may be 
entered at any session, and consecutive courses will follow in successive summers. 
- The courses to be offered in statistics in 1960 at the University of Florida 
are as follows: Statistical Methods I, II and III, through Sample Survey Methods; 
Statistical Theory I, II and II, including Probability, Inference and Least 
Squares; Statistical Problems; Advanced Statistical Inference; and Response 
Surfaces. In addition, a number of courses in the Mathematics Department 
will be available. 

The National Science Foundation is making available to the University of 
Florida grants for college teachers of statistics and college teachers of mathe- 
matics who wish to attend the 1960 session. The stipend is $75 per week for the 
six weeks of the session plus additional amounts for dependents and travel 
allowances. Tuition will also be paid by the National Science Foundation. 
Participants must meet all the admission requirements of the Graduate School 
of the University of Florida and must be admitted thereto or must be a graduate 
student in good standing at one of the cooperating institutions. Applicants for 
these grants should be employed by an institution of higher learning as a teacher 
of mathematics or statistics; those from institutions wherein there is no oppor- 
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tunity for formal training in modern inferential statistics and probability will 
be given priority. No geographical or age limitations will be imposed, though 
preference will be given to the younger age group if other things are equal. 
The previous academic and professional record of applicants will be considered. 
Those who received grants in a previous session who meet the other conditions 
of eligibility are encouraged to apply. Applications for grants should be post- 
marked not later than February 15, 1960 to be assured of full consideration. 
Requests for application blanks for the summer session and for National 
Science Foundation grants should be addressed to Dr. Herbert A. Meyer, 
Statistical Laboratory, University of Florida, Box 3568, Gainesville, Florida. 


GRADUATE ASSISTANTSHIPS IN EXPERIMENTAL STATISTICS 
at NortH Caro.ina State COLLEGE 


Graduate Program 


The Department of Experimental Statistics of North Carolina State College 
offers programs of study leading to the Master of Science and Doctor of Philoso- 
phy degrees. Since 1947, the College has granted 48 Master’s degrees and 53 
Doctorates in Experimental Statistics. The recipients of these degrees are now 
employed by medical schools, statistical departments, agricultural experiment 
stations, industrial and private research organizations and military and other 
governmental agencies. 

At the present time the Department of Experimental Statistics has 13 staff 
members actively participating in the graduate program. Special research 
programs have been developed in quantitative genetics, plant and animal 
sciences, sample surveys, industrial development and engineering and the social 
sciences and economics. Research is oriented to the development of basic 
statistical theory and techniques of design and analysis for surveys and experi- 
ments. Close liasion is maintained with the Department of (Mathematical) 
Statistics at Chapel Hill in order to strengthen the offerings in statistical theory. 
There is also a working arrangement with the Department of Biostatistics in 
the University School of Public Health at Chapel Hill, whereby their graduate 
students can major in Experimental Statistics and minor in the Division of 
Health Affairs. The Department of Experimental Statistics has an IBM 650 
Electronic Data Processing Machine, which provides for extensive statistical 
research in empirical methods as well as the processing of large masses of ex- 
perimental data. 


Graduate Assistantships 


Students who have majored in an applied field and who have a minimum of 
one year of calculus, or students who have majored in statistics or mathematics 
are encouraged to apply for graduate assistantships in Experimental Statistics. 
These assistantships carry stipends of $175 to $200 per month. If a graduate 
assistant has a satisfactory course record, he can complete the requirements 
for a Master’s degree, including a thesis, in 2 years (in less time if he takes 
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courses during the summer). A student with a Master’s degree in Statistics 
can complete the requirements for the doctorate in two years. 

Further information regarding graduate assistantships and the necessary 
application blanks can be obtained from the Dean of the Graduate School or 


the Department of Experimental Statistics, North Carolina State College, 
Raleigh, North Carolina. 


All assistantships will be awarded by April 1, 1960. 


SHort Course oN DEsIGN oF EXPERIMENTS 


To be held June 7 to 17, 1960 at the 
Memorial Center, Purdue University, Lafayette, Indiana 


The course is for statisticians, quality control personnel, engineers, and others 
in industry who are concerned with planning and interpreting the results of 
industrial experiments. 


Topics: 


Review of Analysis of Variance—use of random fixed, and mixed models 
in Analysis of Variance—nested models 


Principles of Experimental Design—balancing, randomization, replication; 
etc. 


Variance Component Analysis 

Randomized Blocks and Latin Squares 

Factorial Experiments—2", 3”, etc. 

Split Plot Design 

Confounding in Factorial Experiments 

Incomplete Block Designs 

Fractional Replications 

Box-Wilson Response Surface Techniques for Optimality 


Course Director: 
Charles R. Hicks 


Short Course Staff: 


Clyde Y. Kramer, professor of statistics, Virginia Polytechnic Institute 

Gayle W. McElrath, professor and head, Industrial Engineering Division, 
University of Minnesota 

Irving W. Burr, professor of mathematics and statistics, Purdue University 


Virgil L. Anderson, director, Statistical and Computing Laboratory, 
Purdue University 


Eligibility: 


It is assumed that those who apply have had previous statistical training 
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including work on tests of hypotheses (or significant differences), linear cor- 
relation, and at least an introduction to analysis of variance. 

Familiarity with a text such as Dixon and Massey’s Introduction to Statistical 
Analysis will be adequate. 


For further information write Professor C. R. Hicks, associate professor of 
mathematics and statistics, Purdue University 
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Errata 


“Factorial Experiments in Life Testing” 
Marvin ZeELEN, TECHNOMETRICS, Vol. 1, No. 3, 1959 


Table 4C, p. 276 should read 


TaBLE 4C 
Summary of Tests of Significance 


x? (critical 
M_ | value at 0.05 level) 


Q 


Hypothesis 


a; 

bj 

Cij 
bjcsj 
aiCi; 
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A Method of Studying Manner of Growth S.C. Pearce 
Analysis of Covariance for a 3 X 4 Triple 

Rectangular Lattice Design (3 associate 

p-.b.i.b.) Bernard S. Pasternack 
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H. A. David 
A Statistical Model for Diagnosing Zygosis by 
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Seymour Geisser 
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On a Test for Order J. B. Chassan 
Partitioning the “Between Slopes” Sum of Squares 

for Forest Growth Data R. M. Cormack 
Tables for W. L. Steven’s “Asymptotic Regression” H. Linhart 


Biometrics is published quarterly. Its objects are to describe and exemplify the use of 
mathematical and statistical methods in biological and related sciences, in a form assimilable 
by experimenters. Members of the American Statistical Association may subscribe through 
the Association at the rate of $4.00 yearly. The annual non-member subscription rate is $7.00. 
Inquiries, orders for back issues and non-member subscriptions should be addressed to BIO- 
METRICS, Department of Statistics, Florida State University, Tallahassee, Florida. 
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souree, volume number and page. References to books should include pub- 
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on a typewriter carefully inserted in ink. Authors are asked to keep in mind the 
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